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because this oxidizer's stored aboard 


Nitrogen Tetroxide provides immediate, hypergolic igni- 
tion with amine fuels. Even after 9 years of storage, 
nitrogen tetroxide is ready for instant use. Combustion 
efficiencies closely approach the theoretical 99%. 
N.O,—Requires no refrigeration—can be stored indefi- 
nitely in missiles at launching sites—Can be used with 
most fuels—including those containing carbon.—Elimi- 
nates rough starts—fast reaction prevents accumulation 


BASIC TO 
AMERICA'S 
PROGRESS 


of propellants in thrust chambers.—Allows throttleable 
control of motors. 

We'll gladly supply technical literature, including: a 
59-page Product Bulletin, and a brochure entitled “Large 
Scale Handling of Nitrogen Tetroxide.” 

For specifications and local offices, see our insert in 
Chemical Materials Catalog, pages 475-482 and in Chemi- 
cal Week Buyers Guide, pages 37-44. 


NITROGEN DIVISION 
Dept. NT14-21-1, 40 Rector Street, New York 6, N.Y. 
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We, at the Research Laboratories, are convinced that it is 
worthwhile. For example, one of the fundamental scientific pro- 
grams being undertaken is a study of the mechanism of solidifica- 
tion of polyphase alloys. 

This study is still in its early stages. However, we now have a 
better understanding of the solidification process. And we have 
demonstrated that the microstructure of certain eutectics can be 
significantly altered by controlled solidification. The contrast is 
illustrated above. 

It is possible that these binary alloys with preferred orientation 
may have unusual physical and mechanical properties. 

Perhaps you are interested in corporate—sponsored studies into 
the fundamental nature of matter. If so, we can offer a research 
environment that seems to stimulate scientists to extra achieve- 
ment. It may be due to the way we encourage an internal cross- 
fertilization of ideas; or because of a unique complex of comple- 
mentary services that free the scientist from routine analysis and 
experimental chores. 


400 Main Street, East Hartford 8, Conn. 
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Opportunities for 
scientists in the phys- 
ics of solids, liquids, 
gases, and plasmas. 
Current studies range 
from the fundamental 
properties of matter 
to the application of 
scientific knowledge 
to promising new 
products. 


For more specific information, 
please write io Mr. W. F. 
Walsh, at the Research Lab- 
oratories. 
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Consistently flight 
performance Project Mercury 


confirms unsurpassed 
THIOKOL solid rocket motors. 


Time after time, NASA’s workhorse, Little ’ 
Joe, has soared into space, checking out 
the workability of materials, propulsion and 
escape systems, and reaction of research 
animals to the environment of space flight. 

Pollux, Recruit, Castor—solid rocket mo- 
a tors from THIOKOL’s Elkton and Redstone 
Divisions—have unfailingly provided the 
thrust and power for Little Joe in its devel- 
«opmental flights. 


ti Sag THIOKOL’s record of propulsion reliability 
he in the spatial program is long and brilliant, 
reaching back to the X-17 which flew suc- 
ittle Joe Nas Carrie is researcn an e- 
way cessfully m 96% of its launches, and to velopment capsule and research animals 
eS 2° earlier research vehicles. to varying altitudes to obtain engineering 
has In NASA’s Little Joe series, THIOKOL and medical ppt to = man 
ie have developed up to 250,000 lbs. thrust, used in these missions are virtually off- 


ene today’s ICBM class. Smaller THIOKOL rock- the-shelf items and are available to other 
ets have been used to free escape capsule 
from booster. 


fe Chemical Corporation 

® BRISTOL, PENNA. 

Plants in: 

TRENTON, N. J.; MOSS POINT, MISS.; DENVILLE, N. J.; ELKTON, MD.; HUNTSVILLE, ALA.; 
MARSHALL, TEXAS; BRIGHAM CITY, UTAH. 


®Registered trademark of the Thiokol Chemical Corporation 
for its liquid polymers, rocket propellants, plasticizers, and other chemical products. 
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Space Propulsion for the future... 


Los Alamos Scientific Laboratory has 
the major responsibility for research, de- 
velopment and testing in the AEC-NASA 
Rover program .. . another of the many 
investigations at Los Alamos into peace- 
time uses of nuclear energy. 


PHOTO: First field test of a KIWI 
nuclear propulsion reactor. 


For employment 
i information write: 
| Personnel Director 


Division 60-61 


LOS ALAMOS NEW MEXICO 
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Scope of ARS JOURNAL 


This Journal is devoted to the advance 
ment of astronautics through the disse nina- 
tion of original papers disclosing new scientific 
knowledge and basic applications of such 
knowledge. The sciences of astronautics are 
understood here to embrace selected aspects 
of jet and rocket propulsion, space!light 
mechanics, high speed aerodynamics, flight 
guidance, space communications, atmosp!icric 
and outer space physics, materials and struc- 
tures, human engineering, overall system 
analysis, and possibly certain other scientific 
areas. The selection of papers to be printed 
will be governed by the pertinence of the t \pic 
to the field of astronautics, by the currer: or 
probable future significance of the reses 
and by the importance of distributing th: in- 
formation to the members of the Society ind 
to the profession at large. 


Information for Authors 


Manuscripts must be as brief as the pr per 
presentation of the ideas will allow.  !cx- 
clusion of dispensable material and con. se- 
ness of expression will influence the Edi: ors’ 
acceptance of a manuscript. In term. of 
standard-size double-spaced typed page a 
typical maximum length is 22 pages of «xt 
(including equations), 1 page of referees, 
1 page of abstract and 12 illustrations. 
Fewer illustrations permit more text, and _ ice 
versa. Greater length will be accept: ble 
only in exceptional cases 

Short manuscripts, not more than one 
quarter of the maximum length stated for ull 
articles, may qualify for publication as 
Technical Notes or Technical Comme ‘ts. 
They may be devoted to new developm: nts 
requiring prompt disclosure or to comm: nts 
on previously published papers. Such ms .u- 
scripts are published within a few month... of 
the date of receipt. 

Sponsored manuscripts are publis ed 
occasionally as an ARS service to the ind is- 
try. A manuscript that does not qualify ‘or 
publication, according to the above-stsied 
requirements as to subject, scope or length, 
but which nevertheless deserves widespread 
distribution among jet propulsion engine: rs, 
may be printed as an extra part of the Journal 
or as a special supplement, if the author or 
his sponsor will reimburse the Soc iety for 
actual publication costs. Estimates are 
available on request. Acknowledgment of 
such financial sponsorship appears as a 
footnote on the first page of the article. 
Publication is prompt since such papers are 
not in the ordinary backlog. 

Manuscripts must be double spaced on one 
side of paper only with wide margins to allow 
for instructions to printer. Include a 100 to 
200 word abstract. State the authors’ posi- 
tions and affiliations in a footnote on the 
first page. Equations and sy mbols may be 
handwritten or typewritten; clarity for the 
printer is essential. Greek letters and unusual 
symbols should be identified in the margin. If 
handwritten, distinguish between capital and 
lower case letters, and indicate subscripts and 
superscripts. References are to be grouped at 
the end of the manuscript and are to be given 
as follows. For journal articles: Authors 
first, then title, journal, volume, year, page 
numbers; for books: Authors first, then title, 
publisher, city, edition and page or chapter 
numbers. Line drawings must be clear and 
sharp to make clear engravings. Use black 
ink on white paper or tracing cloth. Lettering 
should be large enough to be legible after 
reduction. Photographs should be glossy 
prints, not matte or semi-matte. Each illus- 
tration must have a legend; legends should be 
listed in order on a separate sheet. 

Manuscripts must be ac companied by 
written assurance as to security clearance in 
the event the subject matter lies in a classified 
area or if the paper originates under govern- 
ment sponsorship. Full responsibility rests 
with the author. 

Preprints of papers presented at ARS 
meetings are automatically considered for 
publication. 

Submit ripts in duplicate (origi- 
nal plus first carbon, with two sets of 
illustrations) to the Managing Editor, ARS 
i.” 500 Fifth Avenue, New York 


ARS JOURNAL is published monthly by 
the American Rocket Society, Inc. and the 
American Interplanetary Society at 20th «& 
Northampton Sts., Easton, Pa., U. S. 
Editorial offices: 500 Fifth’ Ave. , New York 
36, N. Y. Price: $12.50 per year, $2.00 per 
single copy. Second-class mail privileges 
authorized at Easton, Pa. This publication 
is authorized to be mailed at the special 
rates of postage prescribed by Section 132.122. 
Notice of change of address should be sent 
to the Secretary, ARS, at least 30 days prior 
to publication. Opinions expressed herein 
are the authors’ and do not necessarily refle:t 
the views of the Editors or of the Society. 
© Copyright 1960 by the American Rock«t 
Society, Inc. 
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Multi-Use 
Automated 


The recent demonstration of multi-purpose 
test equipment (MPTE), developed by 
RCA under a series of Army Ordnance con- 
tracts, highlights a new dimension in auto- 
mated multi-use systems support and culmi- 
nates a long-term RCA effort in this field. 
This General Evaluation Equipment is an 
automated, transistorized, dynamic check- 
out system. It contains a completely modu- 
larized array of electronic and mechanical 


Tmk(s) ® 


evaluation equipment, capable of checking 
a variety of electromechanical devices, 
ranging from radar subassemblies to missile 
guidance computers. MPTE provides the 
stimuli, programming, control, measure- 
ment and test functions for the NIKE AJAX, 
NIKE HERCULES, LACROSSE, HAWK 
and CORPORAL missile systems and has 
been extended to other weapons systems 
related to our defense efforts. 
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Prime Contractors to the 
NASA, Air Force, Navy 


Employment Opportunities Available 
for Scientists and Engineers with Above 
Average Abilities 
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When crews of SAC’s Ist Missile Division successfully launched the USAF ICBM Atlas from Vandenberg Air Force Base, September 
9, 1959, the world became aware that the United States had brought into being a formidable retaliatory power for peace. Within four 
months after the first operational launch, the Air Force doubly underlined this missile’s capability. On a single day, January 26, 
1960, the 16th and 17th consecutive successful Atlases were fired intercontinental ranges to predetermined targets from both At- 
lantic and Pacific bases. 

After only five years of intensive development, including concurrent research, testing and fabrication under this nation’s top mil- 
itary priority, Atlas is extremely versatile as well as powerful. It was the Project Score satellite vehicle and is scheduled for use in 
Project Mercury, the Man in Space Program, and in other space exploration missions. Thus, used as a booster for space projects, 
Atlas provides the nation with a key capability in scientific as well as military applications. 

Space Technology Laboratories provides the systems engineering and technical direction for the Atlas as well as other portions 
of the Air Force Ballistic Missile Program. Much of what was learned in building Atlas has helped cut the lead-time in the develop- 
ment of such other Air Force Ballistic Missiles as Thor, Titan and Minuteman. 

Among the industrial organizations which have worked in concert in developing Atlas are such major contractors as: Convair, 
Division of General Dynamics Corp. for airframe, assembly and test; General Electric Co. and Burroughs Corp. for radio guidance; 
Arma, Division of American Bosch and Arma Corp. for inertial guidance; Rocketdyne Division of North American Aviation, Inc., for 
propulsion; General Electric Co. for re-entry vehicle; Acoustica Associates for propellant utilization. 


America’s first 
intercontinental ballistic 
missile...is helping to 
bear the burden of today’s 
power for peace 


The continuing development of Atlas as well as other USAF missiles and related space probes, has created impor- 
tant positions on STL’s technical staff for scientists and engineers with outstanding capabilities in: thermody- 
namics, aerodynamics, electronics, propulsion systems, structures, physics, computer technology, telemetry, and 
instrumentation. If you believe you can contribute in these or related fields and disciplines, you are invited to send 
your resume to: 


SPACE TECHNOLOGY LABORATORIES, INC. > / 


P. 0. Box 95004, Los Angeles 45, California, Attention: Richard A. Holliday Los Angeles * San Diego 
Santa Maria * Sacramento Denver Cheyenne Cape Canaveral Washington, D.C.» Manchester, England Singapore Hawaii 
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ARS Expands Its Publication Program 


T HAS BEEN well known for some time, by 
ARS members and by contributors to the 

ARS Journat, that the flow of manuscripts and 
preprints to the JourNnat has been running far 
above the capacity of the JouRNAL to publish 
them. This problem has been a matter of con- 
cern to the Society’s publication staff and its 
Board of Directors for at least a year. 

The Board has always regarded the publication 
of technical papers as one of the most fundamental 
duties of our Society, and it has always voted 
generous budgets for the JouRNAL, in comparison 
with other engineering societies of comparable 
size. The JourNAL budget programmed for this 
year is larger than that of last year, and new plans 
have been approved for tapping additional 
sources of revenue to permit even greater expan- 
sion of the JouRNAL next year. 

Despite this expansion program for the ARS 
JouRNAL, which will always be limited by the 
revenues available to the Society, a statistical 
forecast indicates that there will still be an excess 
of qualified papers above the publishing capacity 
of the Journat. To solve this problem, a new 
publishing vehicle is needed that can pay its own 
way and expand or contract according to the 
flow of papers. In particular, a severe load is 
expected as a result of the increasing number of 
specialized symposia in various fields that are 
being scheduled by the Technical Committees of 
the Society for 1960, 1961 and future years. 

A solution has been found. The ARS has 
signed a long-term contract with Academic 
Press, Inc., for the publication and distribution of 
a series of books that will function as a supple- 
ment to, and an extension of, the ARS JouRNAL. 
The first book in the series will comprise a 


selection of the papers of the Solid Propellant 
Rocket Research Symposium, and will be placed 
on the market about Sept. 15, 1960. 

Some of the features of this new program are 
the following: 


1 Each book will be cloth-bound, hard- 
covered, good looking, about 63 X 9% in. in 
outer dimensions. Each book will be a dignified 


addition to the owner’s library. 


2 The title of the series will be: “‘ Progress in 
Astronautics and Rocketry.” The subtitle will 
be: ‘‘A Publication of the American Rocket 
Society.” The first volume will be titled: 
“Volume 1: Solid Propellant Rocket Research,” 
and underneath this line will be printed: “* Based 
on a Symposium held at Princeton University, 
January 28-29, 1960. Each book will be edited 
by a specialist (or specialists) in the field, and his 
name (or names) will be featured as editor (or 
(editors). 


3 The price will be down around $5.00 for a 
book about 500 pp. in length, which will help to 
achieve the widest possible distribution. Similar 
low prices will be established for future volumes. 
For the Society, this will be a nonprofit venture. 
Our aim is to serve the technical community. 


4 The Series will not be restricted in concept 
simply to the proceedings of ARS specialized 
symposia. Not all symposia held by the ARS 
will be published in this form. The contents will 
be selective; that is, the Series will be governed 
by criteria of acceptance similar to those of the 
ARS Journat, although subjects of a nonarchive 
character may be covered by individual books 
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from time to time. A very important rule is that 
eac!: volume will be coherent in itself as to subject 
maiter. Heterogeneous collections will not be 
slapped together. The purchaser will not be 
expected to buy a lot of unrelated material in 
order to get the papers he wants. 


5 Promotion of sales will be vigorous, so- 
phisticated, at frequent intervals, and, we believe, 
effective. Our aim is to achieve a 3000 to 5000 
sale of each volume of the series, as a minimum. 
Libraries will be encouraged to place standing 
subscription orders for the entire series, and to 
maintain complete sets right from Volume 1, 
just as they do for periodicals. Low price and 
vigorous selling are important for the achievement 
of broad distribution of the papers to the inter- 
ested scientific community, both in the U. S. 
and abroad. 


6 The projected time lag from the date of a 
meeting to the date of distribution of the printed 
book will be about three months. We expect 
this to apply to the July ARS Symposium on 
Liquid Rockets and Propellants, and to all future 
collections and symposia. Release will occur 
while the subject is still “hot.” Reprints can be 
supplied to the author without delay, at a moder- 
ate price. (The number will be limited in order 
not to undermine the market for the book un- 
duly.) 


These features of the Series, that is, dignified 
book style, sequential publishing, low price, 
prompt release, availability of reprints, critical 
standards of acceptance, adaptability to non- 
proceedings collections, coherent contents, and 
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vigorous, repetitive selling, all combine to make 
this Series attractively different from the run- 
of-the-mill ineffective proceedings publications 
that are sold from time to time by various 
groups. This project has been approved with 
enthusiasm by the Technical Committee Chair- 
men of the ARS. ‘ 

As a result, about five volumes are planned 
for 1960, and still more for 1961, in different 
fields, of course. 

To achieve the low price and prompt publica- 
tion, letterpress printing will have to be avoided. 
Instead, offset printing will be employed, with 
plates photocopied from the authors’ typed sheets. 
For this purpose, all authors will be asked to 
prepare their manuscripts according to special 
typing instructions on photocopy sheets that will 
be supplied by ARS. Although this may seem 
to be an additional burden on the author, we feel 
it is a small price for him to pay for the advantage 
of prompt publication of his paper in a dignified, 
widely advertised book. 

If this plan is successful, we can expect to see 
several results flowing from it. One will be a 
reduction in the ARS Journat backlog, a reduc- 
tion in publication delay, and as a consequence, 
an increase in offerings of worthy manuscripts to 
the Society. 

A second result will be greater readiness on the 
part of our Technical Committees to hold special- 
ized symposia, since they were often deterred in 
the past by the insoluble problem of getting the 
papers of such symposia published. Finally, 
with considerable optimism, we expect that this 
program will prove so successful that other tech- 
nical societies will adopt similar plans to solve 
their publication backlogs. 


Martin Summerfield 
Editor, ARS JOURNAL 
Editor, ARS Progress Series 


Contract Signed with Academic Press for ARS Progress Series 
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Radio Astronomy 
Observations From Space 


FRED T. HADDOCK' 


University of Michigan 
Ann Arbor, Mich. 


A description is given of four radio astronomy experiments that can be carried out with simple 
very low gain antennas at frequencies which are too low to penetrate, free of distortion, Earth’s 
ionosphere. These are: 1 Measurement of the integrated radio flux from the galaxy at frequencies 
below 10 mc; 2 measurement of the dynamic spectra of solar bursts of frequencies below 30 me; 
3 measurement of the dynamic spectra of bursts at frequencies below 30 mc from the planet 
Jupiter; 4 measurement of the distribution of electron densities outward from the F-layer of 
Earth’s ionosphere. These experiments can best be conducted aboard a highly eccentric artificial 
Earth satellite with an apogee greater than 5 to 10 Earth radii. 


URING the past decade or so extensive radio astronomy 
observations have been conducted at many frequencies 
ranging principally from 10 to 10,000 me ( a frequency of 10 
me corresponds to a vacuum wave length of 30m). Earth’s 
lower atmosphere strongly absorbs and emits radio waves 
above a frequency of 10,000 me, thereby limiting higher 
frequency observations. But since this observational limita- 
tion can be largely overcome by flying the receiving equipment 
in balloons or high altitude aircraft, satellite experiments at 
these frequencies are not yet urgent. At frequencies from 
about 10 to 10,000 me, the vertical atmosphere is normally 
highly transparent. However, the ionosphere and troposphere, 
by distorting the phase fronts of incoming waves, limit the 
measurement of size and position of radio sources and the 
radio emission detail that can be seen on the sun and planets. 
Satellite observations would also overcome this limitation in 
resolution, but would require very large antennas or inter- 
ferometers; one would need an antenna system extending 
over 10 km to resolve 1 min of are at 100 me. 

On the other hand, satellite observations at frequencies 
below about 30 me, even with low gain antennas, would be 
extremely valuable. It would be possible to extend the 
dynamic radio spectra of solar bursts and bursts from Jupiter 
down to 5 mc or lower and to check and extend the work of 
Reber and Ellis (1)? who extended the spectrum of the general 
cosmic radio background emission to 1 me by observing 
through occasional “holes” in the ionosphere. It would also 
be possible to determine the free-electron density distribu- 
tion in Earth’s atmosphere above the F-layer of the iono- 
sphere and determine how this distribution blends into that 
of interplanetary space or the solar corona. These are im- 
portant and basic experiments which hold promise of greatly 
increasing our knowledge of solar flares, solar corpuscular 
streams, the sun’s outer atmosphere, the atmosphere of 
Jupiter, the interstellar medium, the galactic corona and the 
origin of cosmic rays, as well as of providing valuable design 
data for future satellite experiments relating to ambient radio 
noise levels (both steady-state and transient), the characteris- 
tics of the radio transmission through the solar system, and 
the leakage of natural static and manmade signals up through 
the ionosphere. 


Presented at the ARS Semi-Annual Meeting, June 8-11, 1959, 


San Diego, Calif. 
1 Professor of Astronomy and of Electrical Engineering. 
? Numbers in parentheses indicate References at end of paper. 
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When more elaborate and highly directional satellite antcn- 
nas are developed, it will be possible to extend the radio 
spectrum of the brighter radio sources. The determination 
of the spectra of radio sources at low frequencies is required 
to specify the physical parameters in the radiation process 
in the variety of interesting radio sources, such as the colliding 
galaxies in the constellation of Cygnus, the peculiar nebulosity 
in Cassiopeia, the supernova remnant in Taurus, the galactic 
center source, etc. In short, satellite observations make it 
possible to extend the radio astronomy spectrum downward 
by several octaves, to a limit determined only by the free 
electron in the interplanetary medium or solar corona. 

The following radio astronomy satellite experiments can 
be made with simple long wire, short stub, or loop antennas 
with conventional radio receivers and will now be discussed 
individually: 

1 Dynamic spectra of solar radio bursts. 

2 Dynamic spectra of radio bursts from Jupiter. 

3 The spectrum of cosmic radio background radiation. 

4 The distribution of electron densities outward from the 
F-layer. 


Dynamic Spectra of Solar Radio Bursts 


A description of solar activity and flares, how they are 
related to events on the sun and Earth, and their scientific 
and practical importance appears in many texts (2). Certain 
of the greater solar flares are accompanied by large bursts of 
radio emission which begin at the higher frequencies 
(above about 200 mc) and appear progressively later at the 
lower frequencies. The observed rates of this downward 
drift in frequency fall into two broad but distinct classes, as 
first noted by Wild and McCready (3) in 1949 from their 
sweep-frequency measurements of solar radio bursts. They 
denoted the fast drift rate bursts [of the order of 20 (mc/s) 
sec~!] as type III and the slower drift rate bursts [of the 
order of 0.25 (me/s) sec] as type II bursts. They denoted 
as type I the narrow spectrum (a few me wide), short duration 
(less than 1 to 20 sec) impulsive “noise storm” bursts. 
Wild (4) has interpreted the type III and II bursts as caused 
by corpuscular streams shot out from the region of a solar 
flare, at velocities of the order of 50,000 km per sec for type 


III and 1000 km per sec for type II. The streams excite 


plasma oscillations of decreasing frequency as they pass 
through successively more rarefied layers of the corona. 


ARS JourNAL 


cosm: 
highe 
event 
Th 
radio 
ionos| 
xX 10) 
mum 
F-laye 
perige 
have | 
cessiol 
for me 
Fig. 
at Ea 
tensiti 
turbed 
siopeia 
backgr 
It ca 
greater 
gain at 
log-per 
directe 
axes, W 
thereby 
Howev 
stand t 
ment a 


1 


TI 
| for 
he 
| ma 
vel 
sto 
ty] 
cor 
on 
der 
a 
| (be 
tur! 
typ 
low 
cre: 
| rela 
im} 
belc 
whi 
| sph: 
F reqi 
tion 
Jupi 
| mak 
for t 
sun, 
| the | 
dista 
7 dece 
its fi 
the | 
the b 
the 1 
secor 
It. 
= will | 
corre 
& 


The approximate heights of different plasma levels are known 
for optical eclipse determinations of coronal electron densities; 
hence the outward radial component of velocity can be esti- 
mated from the frequency drift. Wild further noted that the 
velocities of particles presumed to cause great geomagnetic 
storms following great flares are of the same order as the 
type II burst velocities, and that the time of flight velocities 
computed for the rare solar cosmic ray increases following 
very great flares correspond to the type III burst velocities. 

Not all radio bursts of type II or III give rise to a 
magnetic storm or a cosmic ray increase. There is some evi- 
dence that a number of type II and type III events start at 
a high frequency, reach only a certain limiting frequency 
(because of a lack of kinetic energy in the stream?) and are 
turned back or stopped. It has been tentatively noted that 
type II and type III bursts which appear to terminate at a 
low frequency are not associated with geomagnetic index in- 
ere:ses (5). In general, the type II bursts do show some cor- 
relation with the geomagnetic increases. Therefore it is 
important to make observations of solar spectra in a band 
below about 20 me. This is difficult, even at those frequencies 
which normally penetrate the ionosphere, because the iono- 
sphere is disturbed at the very time these observations are 
required. The satellite data can be correlated with observa- 
tions from the ground to help identify the particular event 
observed and also to help identify bursts from the planet 
Jupiter and terrestrial interference. 

The importance of this experiment is twofold: It would 
make possible the further checking of the stated hypothesis 
for the cause of the great bursts of radio emission from the 
sun, and it would also make possible, under this hypothesis, 
the determination of the decrease of electron density with 
distance from the sun, and perhaps the acceleration or 
deceleration of the corpuscular stream in the initial phase of 
its flight from the sun to Earth. 

The observation could consist, for example, of measuring 
the burst intensity as a function of frequency and time over 
the band of 5 to 30 me with a receiver band pass of 0.2 me while 
the receiver is tuned over the band one or more times a 
second. 

It is expected that some type III and type II burst spectra 
will be recorded over the entire 5- to 30-mc band and would 
correlate with geomagnetic storms, aurorae and perhaps solar 
cosmic ray increases. Some will be detected only at the 
higher frequencies and may not be correlated with terrestrial 
events. 

The orbit perigee should be sufficiently high to eliminate 
radio propagation disturbances of the solar signal by Earth’s 
ionosphere. The critical electron density for 5 me is 3.1 
x 10°em~%. This is about 10 times smaller than the maxi- 
mum over England at noon during November 1957. ‘The 
F-layer is likely to interfere with observations unless the 
perigee is above about 1200 km. A desirable orbit would 
have a polar inclination which would result in an orbit pre- 
cession that would keep the sun under continuous surveillance 
for months. 

Fig. 1 displays the variation with frequency of flux density 
at Earth (in MKS units) for average peak solar burst in- 
tensities (5,6). The flux density from the quiet or undis- 
turbed sun, from the strongest cosmic radio source (Cas- 
siopeia-A), and from the integrated intensity of the cosmic 
background radiation (labeled G=1) is also shown (1). 

It can be seen that the average burst level is about 10 times 
greater than the cosmic background noise level when a low 
gain antenna is used. A long wire Vee antenna, perhaps of 
log-periodic structure, having its beam fixed in space and 
directed toward the sun from a satellite stabilized about three 
axes, would give a large increase in signal to noise ratio, and 
thereby make is possible to detect weaker solar radio bursts. 
However, such an antenna must be sufficiently stiff to with- 
stand the tidal force of Earth in order to prevent its align- 
ment along the Earth-satellite radius vector, unless the 
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Fig. 1 Radio flux density spectrum of average solar burst in- 

tensities extrapolated downward from about 30 mc, the galactic 

background flux density measured on a spinning dipole pattern 

and the spectra of the strongest radio source Cas A and the quiet 
sun 


satellite is in a planetary orbit around the sun. The maximun 
acceleration due to the tidal force near Earth is about 5 X 
10-8 1 of surface gravity, where | is the distance (in meters) 
from the satellite center of gravity. 

It would be important to have a number of related ex- 
periments operating concurrently in the same satellite and on 
the ground. For example, experiments for measuring UV, 
x-rays, cosmic rays, geomagnetic stream particles, Lyman-a 
spectroheliograms, geomagnetic perturbations, etc., could be 
in the same vehicle, or in other satellites, while the standard 
IGY experiments are conducted more intensively on the 
ground. 


Dynamic Spectra of Radio Bursts From Jupiter 


In 1955, Burke and Franklin (7) unexpectedly discovered 
emission of intense bursts from the planet Jupiter at 22 mc. 
Since then a number of very interesting and important facts 
have emerged (8). It has been known for several years that 
about 20 per cent of the time the planet Jupiter emits very 
intense bursts of radio emission in the frequency band of 14 
to 27 me. Each individual burst has a rather narrow band- 
width of the order of 1 me. They occur independently at 
different frequencies and last for 0.1 to 15 sec. The bursts 
are most numerous at about 18 mc; observations below 18 
me are rare, perhaps because of interfering signals and dis- 
turbance or blocking by Earth’s ionosphere. The detec- 
tion of bursts at frequencies above 30 mc has not been con- 
firmed. Therefore the spectra of the radio bursts from Jupiter 
have an extremely sharp high frequency cutoff. Nothing 
definite is known about the low frequency cutoff; however, 
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if a Jovian ionosphere controls the duration of burst activity, 
then low frequency cutoffs from 8.5 to 20 mc are indicated. 
On the other hand, this concept has been criticized by Carr 
et al. (8). Gallet has suggested that corpuscular streams are 
ejected from the “surface” of Jupiter through a Jovian iono- 
sphere, creating plasma oscillations and radio bursts (8). 

A fact bearing on the importance of the radio bursts from 
Jupiter is that the radio energy of a single large burst is 
about 10'? watt-sec or of the order of the energy expended 
in an average volcano. The overall significance of radio 
bursts from Jupiter can be appreciated when it is noted that 
the radio observations to date strongly suggest the existence 
of an ionosphere on Jupiter, a magnetic field with a strength 
of at least 5 gauss, or 10 times that of Earth, and that the 
energy sources for the radio bursts are fixed to the surface of 
a solid rotating body. Recently it has been shown by Carr 
et al. (8) that the period of rotation of this body is 11.8 sec 
shorter than the period observed in the nonequatorial atmos- 
phere of Jupiter. The determination of the energy source 
of bursts from Jupiter may bear on the composition and physi- 
cal state of the invisible solid body of Jupiter, and detailed 
dynamic spectra studies may disclose information on the outer 
atmosphere, ionosphere and magnetic field of this planet. 

The observation of Jupiter could employ the same receiver 
and antenna as proposed for the solar burst experiment. 
It is expected that numerous bursts will be found over the 
band of 18 to 30 me. However, knowledge of the range and 
frequency of occurrence of bursts below 18 me will provide 
new and valuable data for testing hypotheses of the emission 
process and influence of a Jovian ionosphere on the bursts. 
Correlation of this range and frequency with solar induced 
effects in Earth’s ionosphere can also be studied. It will 
be especially interesting to observe Jupiter during a sudden 
disturbance of Earth’s ionosphere. 
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Fig. 2 Radio flux density spectrum of Jovian bursts, intensities 

extrapolated downward from about 50 mc, the galactic back- 

ground flux density measured on a spinning dipole pattern and 

the spectra of the strongest radio source Cas A and the quiet 
sun 
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Fig. 2 is the same as Fig. 1 except that the solar burst 
spectrum has been replaced by the Jovian burst spectrum. 
Note that the intensity of Jovian bursts is somewhat less 
than that for solar bursts. Increased antenna gain would be 
desirable, but the orientation of a directive antenna toward 
Jupiter adds complexity. 


Spectrum of Cosmic Radio Background Radiation 


The determination of the low frequency radio frequency 
spectrum of the cosmic background radiation at low and 
high galactic latitudes is required to separate the compone: ts 
of radio emission (both galactic and extragalactic), and 
thermal and nonthermal emission. The thermal and non- 
thermal emission processes that contribute to the cosniic 
background spectrum should produce a maximum brig! t- 
ness temperature at a frequency in the band of 1 to 20 nc; 
the spectrum then rapidly decreases at lower frequencies (‘'). 
Because of the intervention of Earth’s ionosphere, it is not 
possible to make reliable observations from the ground beir- 
ing on this critical point; observations made from a satellite 
would make possible direct observation without absorption 
and reflection by Earth’s ionosphere. The principal valve 
of the experiment is the determination of the important :- 
rameters involved in the process of cosmic background emis- 
sion. The source of the cosmic background radiation is i1n- 
portant from the point of view of cosmolgy, galactic struc- 
ture, the origin of radio sources and cosmic rays. 

The only information available on the cosmic backgrouid 
intensity below 9 mc is from observations made by Reber aid 
Ellis (1) in Tasmania at times when the ionosphere was trans- 
parent (at times of low solar activity and at night). They 
were able to obtain estimates of the brightness of the cosmic 
background in the range of from 1 to 10 me. However, thie 
accuracy of their measurements depends upon estimating the 
size of the transparent hole in the ionosphere. Therefore 
it is proposed that experiments with eight single-frequency 
receivers operating on separate or the same antenna be per- 
formed. The frequencies should cover uniformly the range 
of about 65 ke to about 6 me. This would make possible 
the determination of electron densities over a wide range of 
values and also make possible a great stride forward toward 
extending information on electron densities in the solar corona 
at large distances from the sun or in the interplanetary 
medium. The observation would consist of the continuous 
recording of the intensity on each channel. It would be 
desirable to scan the antenna beam over a great circle per- 
pendicular to the galactic plane, and to continue the observa- 
tions to note changes at the time of solar events and to record 
bursts from the sun and Jupiter. 

We would expect of find that the intensity of the cosmic 
background emission peaks at higher frequencies in the ga- 
lactic plane than at the poles because of absorption by inter- 
stellar electrons which are concentrated in the galactic plane. 
The peak intensity of emission will occur near 0.5 me in the 
plane, depending on the antenna directivity; whereas at the 
galactic pole the peak of emission may occur between 0.02 
and 0.2 me. 

It appears likely that interplanetary electrons will prevent 
an unobstructed view of the galaxy at the lowest frequencies, 
even with the orbit at heights of a few Earth radii. Ifa 
density of 600 em~* exists around the Earth-moon system, 
then the galaxy is blocked from view below about 0.22 me. 
An orbit inclination perpendicular to the plane of the Milky 
Way is desirable only for a directive beam measurement. 
Some directivity is obtained by Earth occulting part of the 
galactic background. A low directivity measurement of the 
radio background intensity can be accomplished with a small 
loop or stub antenna, since antenna radiation efficiency is 
not too important because of the high background noise 
level. 

A directive beam measurement of background radiation 
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can be made by using a very long wire antenna (at least 10 
wave lengths in each half) released from a satellite stabilized 
in three axes, orientated and locked onto Earth, in an orbit 


more or less perpendicular to the plane of the Milky Way—the 
two halves of the antenna being aligned by tidal forces with 
the Earth-satellite vector, one toward Earth and one away, 
each 45-km long (90 lb of AWG No. 30). This antenna 
would have a directive gain of 9 at 65 kc and somewhat higher 
at the higher frequencies. It is possible to obtain a sufficiently 
low ohmic-to-radiation-resistance ratio if a receiver noise 
figure is less than 3 db. The natural “pendulum” frequency 
of this long wire antenna is independent of length and is 
sufficient to keep the wire aligned with Earth and the satellite. 
A two-wire, balanced-fed antenna is required because the 
electrical capacity of the satellite is too low at act as a ground 
plane for efficient operation. 

lig. 3 displays the spectrum of cosmic radio background 
intensity, in temperature units, as extrapolated from higher 
freyuency measurements (10) and as measured at a number of 
fre,uencies from 1 to 4 me by Reber and Ellis (1). The 
sha led area represents speculation of a decrease in intensity 
at low frequencies. The decrease can be intrinsic to the 
emission process or, more likely, to absorption or reflection 
by the electrons in the sun’s outer corona, by interplane- 
tary electrons, or an extension of Earth’s ionosphere. 

It is expected that a sudden drop will be found in the 
background brightness around 0.14 to 0.3 me. 


Distribution of Electron Densities Outward From 
the F-Layer 


lor the radio astronomer, knowledge of the electron 
density surrounding Earth is important, since electrons pre- 
vent cosmic radio waves from reaching an Earth satellite at 
all frequencies below a certain critical value. If an estimate 
of 600 electrons per cm! is correct, it is not possible to observe 
the galaxy, or even the planets, with frequencies below 220 
ke. 

The electron density distribution around Earth is impor- 
tant for a number of geophysical and astrophysical problems. 
It bears directly on studies of Earth’s ionosphere and upper 
atmosphere. It is important for studies of Earth’s magnetic 
field and the interplanetary magnetic field, for the nature, 
origin, and heating of the sun’s corona, the conduction of heat 
from the sun to the interplanetary medium, and the nature 
of zodiacal light. 

The electron density distribution above the F-layer can be 
estimated from ionospheric theory, from the Faraday rota- 
tion of radio waves from moon radar echoes and from radio 
transmissions by rockets or satellites. These determinations, 
however, involve, or are affected by, the dense regions of 
ionozation in the F-layer, and therefore are not as attractive 
as a direct determination of the density surrounding a 
satellite as it travels from the dense to the distant tenuous 
regions. 

A Langmuir probe method for the determination of elec- 
tron temperature and density has been demonstrated on 
rockets in the E-layer of the ionosphere and has been consid- 
ered for Earth-satellite measurments (11). This method 
should be compared with the following suggestion for meas- 
uring ambient electron densities; this is an r-f probe technique 
which may have certain advantages over the d-c Langmuir 
probe. 

This experiment consists of tuning a radio frequency 
oscillator, attached to a matched long wire antenna extending 
from the satellite, over a frequency band from 20 ke to about 
10 me at a rate determined by the rate of change of altitude 
of the satellite. For a circular orbit this sweep rate could be 
rather slow; for an eccentric orbit the sweep period should be 
short compared to the time it takes a satellite to change its 
altitude. As the oscillator is swept over the frequency range, 
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Fig. 3 The integrated cosmic ray background spectrum extrap- 

olated below 10 mc with a few points measured by Reber and 

Ellis. Shaded area indicates probable spectrum at lower 
frequencies 


the change in the impedance of the oscillator circuit is re- 
corded. When the oscillator frequency approaches the 
critical frequency of the ambient electrons, there would be a 
change in the record. With magnetic fields at oblique 
angles with respect to the antenna, analysis would be needed 
to interpret the results. It may be possible that changes in 
aspect of the antenna with Earth’s field will give additional 
information about the field. The rate at which information 
must be telemetered to Earth for this experiment depends 
largely on the orbit. The only data needed would be meas- 
urements of the change in the impedance with position of 
the satellite, and perhaps the antenna orientation. It may 
also be useful to record the phase of impedance as well as 
its magnitude. It would be extremely important to deter- 
mine the change in magnitude and distribution of the electron 
density, following a flare or radio burst and during a geo- 
magnetic storm or aurora. 

It is expected that the electron density will be found to 
decrease from the F-layer maximum near 10° cm~* (de- 
pending on local time, season, solar activity and latitude) to 
several thousand cm~ at an altitude of 1000 to 4000 km, and 
then gradually to decrease to a few tens or hundreds at 
several Earth radii. The perigee should be about 250 to 
300 km and the apogee several Earth radii. A polar orbit is 
preferred in order to determine the effect of the geomagnetic 
field on the electron distribution. 

Fig. 4 illustrates schematically some possible distributions 
for the electron density about the F-layer. The right-hand 
ordinates are the critical plasma frequencies for the corre- 
sponding electron densities on the left-hand ordinates. The 
dashed line on the left represents the decrease above 
the F-layer maximum. Its slope is based on an article 
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on Sputnik I and II (12) and a recent Russian paper (13). 
The solid curve is based on a model by Dungey (14) which 
represents the electron density decrease with the distance 
R from Earth’s center by 600 exp (2.5 R./R) cm~* where 
R, is the radius of Earth. The two levels on the right-hand 
side of Fig. 4 labeled H represent the range of values given 
by van de Hulst (15) from a study of the polarization of 
zodiacal light. Six hundred cm~ is the asymptotic value 
by Dungey’s model and was adjusted to this in order to fit 
both the zodiacal light determination and the results obtained 
from magneto-ionic-duct propagation studies (‘“‘whistler’’) (16). 
Later studies of ‘whistlers’’ indicated an electron density of 
300 cm~? at a distance of 3 Earth radii (16). 

Two recent models of the electron distribution of the 
terrestrial “corona” also based in part on very low frequency 
emission data are shown, one by Gallet (17) and one by 
Johnson (18). 

The divergency of these curves reflects our lack of know]l- 
edge which should soon be reduced by satellite experiments. 
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Simulation of Fuel Sloshing 


Characteristics in Missile Tanks 


by Use of Small Models 


H. NORMAN ABRAMSON? 
GUIDO E. RANSLEBEN Jr.* 


Southwest Research Institute 
San Antonio, Texas 


Similitude theory is applied to the problem of fuel sloshing in accelerated tanks to establish 
criteria for the design of model experiments. It is found that dynamic modeling is possible even if 
liquid viscosity is considered, and the ranges of significant parameters and the selection of model 
liquids are discussed. The results of experiments made with small models are compared with those 
obtained with full-scale tanks, for two different types of damping devices. 


LOSHING motions of fuel in partially filled tanks may be 

established by essentially lateral accelerations of the ve- 
hicle. Such fuel motions are of importance, for there is the 
possibility of extreme oscillations if the excitation frequency 
is in the neighborhood of one of the natural frequencies of the 
liquid fuel. The forces exerted by the fuel on the tank walls 
may then lead to perturbations in the flight trajectory, or may 
impose severely high stresses on structural components. 
There is further the ever-present danger of resonance at the 
control system frequencies. 

Suppression of sloshing modes has been attempted by means 
of various mechanical devices and baffle systems, and, there- 
fore, it is of some importance to evaluate the effectiveness of 
such devices. Purely mathematical approaches to this prob- 
lem are exceedingly difficult because of the boundary condi- 
tions introduced by such suppression devices, and are re- 
stricted, therefore, essentially to simple cases. On the other 
hand, it may be possible to evaluate suppression devices by 
means of suitably designed and conducted model tests. 

The problem considered here is the specification of criteria 
for model tests which will give useful quantitative information 
concerning sloshing and slosh suppressors in full-scale fuel 
tanks. This is accomplished by the application of similitude 
theory to obtain appropriate expressions for the model-proto- 
type relationships. These may be employed in turn to yield 
a rational design of model tests. This is in contrast to pre- 
vious experimental studies of the fuel sloshing problem (1-4),4 
which have apparently been conducted without regard to 
similitude requirements, and consequently are believed to be 
of limited usefulness for evaluation of the effectiveness of slosh 
suppression devices. 

The purposes of the present paper are to discuss the basic 
requirements for model studies of fuel sloshing in rigid tanks 
by application of similitude theory, and to show correlation 
between tests made on small models and on full-scale tanks, 
for two different suppression devices. 
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Similitude Relations for Fuel Sloshing 


General 


The vehicle is considered to be in vertical flight with con- 
stant acceleration directed along the flight path, and sloshing 
is assumed to be excited by time-dependent accelerations 
acting essentially normal to the flight path. For purposes of 
the present study, the exciting accelerations will be con- 
sidered either as purely translational or purely rotational 
about an axis normal to the flight path. 

The essential features of the sloshing motion are further 
governed by the physical properties of the liquid and the 
pressurizing gas. The importance of the liquid density and 
total mass is obvious; however, viscosity and surface ten- 
sion require additional comment. If suppression devices are 
incorporated, it is clear that their effectiveness must, to a 
very substantial degree, depend upon the damping forces they 
can provide as the result of viscous action of the liquid. If 
models of very small size are employed, it is conceivable that 
the effects of sloshing may be substantially modified by sur- 
face tension forces acting on the suppression devices and at the 
tank walls. As will be discussed later, it does not seem possi- 
ble to provide rigorous modeling if surface tension forces are 
included, and, further, it is believed that the total contribu- 
tion of such forces would be small compared to the inertial and 
viscous forces. Consequently, surface tension forces will be 
omitted in the analysis. Similarly, it is supposed that the 
forces produced by the pressurizing gas are negligible com- 
pared with the other forces, and therefore will be omitted. 
The gas pressure may be accounted for implicitly, however, 
because it is included in the total pressure at the tank wall by 
direct superposition, since the gas volume remains essentially 
constant during sloshing. 

Further restrictions to the analysis are provided by the as- 
sumptions of small excitations, rigid tank and the existence 
of geometrical similarity in all respects between model and 
prototype. The model and prototype suppression devices 
themselves are considered geometrically similar in all respects, 
and the corresponding inertial and apparent mass forces as- 
sociated with them are neglected in comparison with the vis- 
cous forces. 
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Similitude Relations—Translational Excitation 


The relevant parameters, based on the preceding discussion, 
will be taken as 


= longitudinal acceleration on tank, LT~* 
= tank diameter, L 

= resultant liquid force on tank wall, MLT~? 
= depth of liquid in tank, L 

» = excitation amplitude, L 

= liquid viscosity, 

= liquid density, ML~* 

= excitation period, T 


NOR MP WAS 


An equation relating these eight parameters analytically can 
be written in the general form (5) 


¢'(m™, Te, Tn — m) = 0 {1] 
where 
nm = number of physical parameters involved 
m = number of fundamental dimensions (in this case 


n—m=8-—3=5) 


The z’s are dimensionless combinations of the parameters 
listed. Proceeding in conventional fashion, similitude theory 
yields expressions for the various dimensionless groups. The 
general solution of the form of Equation [1] may then be writ- 
ten in terms of the 7’s as 


Ez h Xo, | 
d dp 


We note that all linear dimensions are scaled in the same ratio 
as the diameter, that the group ar*/d is equivalent to Froude’s 
number, and that the group pd(d/r)/u is equivalent to 
Reynolds’ number. 

Should interest reside not in the resultant liquid force on 
the tank wall, but rather in the resultant liquid pressure, a 
similar analysis would yield 


p = f(F/d*) [3] 


so that in place of the dimensionless force F'/[(pd*)(d/r?) | in 
Equation [2], we may write 


[2] 


pre 
pd® p(d/r)* 


which is in the form of a pressure coefficient (Euler number). 

The preceding analysis has considered the excitation to be a 
displacement, in anticipation of actual model experiments. 
The excitation could, however, have been taken equally well 
as a force P (or even an acceleration), in which case the ap- 
propriate dimensionless group would become 


P/((pd*)(d/r*)| [5] 


Consideration of the surface tension in the foregoing analysis 
would have resulted in a new dimensionless group of the form 
pd*/or*, which corresponds to a Weber number. As will be 
seen later, the simultaneous solution of this new group and 
pd(d/r)/ would be exceedingly difficult, in view of the model 
liquids readily available. Further, we assume that surface 
tension forces are small compared with the inertial and viscous 
forces, and consequently we shall omit any additional con- 
sideration of surface tension. 


[4] 


Similitude Relations—Rotational Excitation 


In this case the pertinent parameters remain the same as in 
the translational case, except that the excitation is now de- 
fined by two parameters, the angular rotation 65 and the loca- 
tion of the rotational axis b. Since we consider only small ex- 
citation amplitudes, it may be noted that b0) ~ Xo, and, 
therefore, the general functional Equation [2] applies to both 
kinds of excitation. The use of a moment M = Pb as the 
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excitation would result in a group formed by Equation [5] 
multiplied by the ratio b/d. 


Modeling Considerations Based on the Similitude 
Relations 


General 


Equation [2], and the associated alternate expressions ‘or 
special cases, governs the design of models to simulate a given 
prototype. If all dimensionless groups in this equation have 
the same numerical value for both model and prototype, thin 
the forces and pressures measured on the model are directly 
applicable to the prototype. 

We note that fluid depth and excitation displacement a:i- 
plitude are scaled in the same proportion as the geometrical 
scale. Thus (the ratio of model to prototype parameters is 
denoted by the subscript r) 


d, = dn/dp 
h, = Xo, = by = d, 0, =1 {7] 


The remaining model parameters are to be established |v 
further study of Equation [2]. 


Modeling Considerations—Viscosity Neglected 


In the case that viscous effects are neglected throughout, 
the Reynolds number group drops out of the analysis given 
previously, and the relation between model and prototype d:- 
pends on the geometrical factors discussed and the time scale 
factor 


= d,/a, 


This result means that any size model in a 1-G acceleration 
field can be used to simulate any size prototype in any ac- 
celeration field by varying the time scale, and would further 
permit the use of any model liquid, since only the density ap- 
pears in the force (or pressure) parameter. 

In presenting measured data, the parameter ar?/d may be 
used as the primary independent variable. Since frequency is 
a more commonly used variable, it may be introduced for the 
inverse of the period to give dw?/a. The force (or pressure) 
group may also be modified by introduction of Equation [8] to 
give 

F Pp 
or 
pad 
The measured data may then be presented with dw?/a as the 
primary independent variable, F/pad* as the dependent 
variable, and the remaining quantities as secondary inde- 
pendent variables. 


[9] 


Modeling Considerations—Viscosity Included 


With viscous forces included, the data may still be pre- 
sented as indicated; however, the model design is now subject 
to certain restrictions arising from the Reynolds number 
group. Introducing Equation [8] we find 


pd(d/r) _ 


[10] 


from which 
d, = (pr/pr)*a, (11) 


This important result shows that the geometrical scale is de- 
termined by the acceleration, liquid density and liquid vis- 
cosity ratios. If it is desired to test models so as to cover : 
range of prototype accelerations, we may either fix the mode! 
diameter and vary the liquid properties, or fix the mode! 
liquid and vary the model size. Since the range of liquid 
properties is very restricted, a combination of these courses 
may be chosen. 
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In order to employ models of reasonably small size to study 
large prototypes, the model liquids must have high density 
or low viscosity, or both, compared to the prototype liquid, as 
inspection of Equation [11] indicates. Mercury, of course, 
has a very high density, but is undesirable for a variety of 
reasons, such as toxicity and cost; some of the organic sol- 
vents have very low viscosity and, therefore, offer some 


promise. 
As an example of model simulation, assume the following 
typical case: 
Prototype Model (methyl- 
Parameter (kerosene) ene chloride) 
Specific gravity 0.83 1.336 
Viscosity 2.5 ep 0.465 cp 
Acceleration 1-5 G 1G 
Diameter 9 ft (to be derived) 


lor a prototype acceleration of 1 G 


For any other prototype acceleration, the model size is 2.14 
divided by the cube root of the acceleration ratio. Thus 


(dm)sq = 2.14/(1/5)'/? = 3.66 ft 


Therefore, a 9-ft diameter prototype with kerosene fuel under- 
going accelerations from 1-5 G may be simulated by a series of 
models of 2.14 to 3.66-ft diameters with methylene chloride as 
the model liquid. The range of excitation frequencies for the 
model is determined from Equation [8] as 

1 

Oo =—- =| = 2.502-0.701 

Tr 

It is clear from the foregoing that the inclusion of an addi- 
tional group to account for surface tension would have 
rendered the model design virtually impossible, because of the 
inability to find a model liquid of appropriate properties to 
satisfy all dimensionless groups simultaneously in a small 
model. 


Modeling Considerations—For Comparison With Full- 
Scale ABMA Tests 


The availability of test results obtained by the Army Ballis- 
tic Missile Agency on a full-scale missile tank provides the 
means for verifying the results of the foregoing analysis. The 
full-scale tank tests were conducted on the ground with water as 
the “fuel’”’ and with the excitation as translational. Two dif- 
ferent types of suppression devices employed in these tests 
were simulated with the model. 

The prototype (full-scale) parameters pertinent to the 
similitude requirements are 


d, tank diameter = 8.75 ft 

p, liquid density = 1.938 slugs/ft* 
M, liquid viscosity = 1.00 ep 

a, acceleration 1G 


The model is necessarily in a 1-G acceleration field, there- 
fore its diameter is determined by Equation [11] as 


m 
pm 1.00 1 


In order to obtain a reasonably small model (4-,'5 scale), 
it becomes necessary to choose a liquid with suitable values 
of viscosity and density, in accordance with this equation. 
Methylene chloride appears suitable as a model liquid, since 
it is nonflammable and relatively free of toxic hazards; how- 
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Z-RING | STIFFENERS 


2 


Fig. 1 Basic tank configuration 


ever, it does have the disadvantage of being a solvent for most 
plastics and thus eliminates the use of models made of trans- 
parent Plexiglas or lucite. Inserting the properties of 
methylene chloride (6) into the previous equation yields a 
model diameter of 4.33 ft, which is somewhat larger than de- 
sirable. Since no other suitable model fluids could be found, 
models of 14.4- and 25.0-in. diameters were finally selected. 
This meant that exact simulation of the ABMA full-scale tests 
(with water) could not be accomplished (although “in-flight”’ 
simulation was possible, as mentioned previously). However, 
it appeared that by employing various model sizes and liquids 
(e.g., water, methylene chloride, methylene bromide), it 
would be possible to provide a sufficient variation in Reynolds 
number to yield useful information regarding the damping 
effectiveness of slosh suppression devices. This will be dis- 
cussed later. 


Experimental Arrangements 


Models 


The models were fabricated from rolled steel cylinders with 
spun conical bottoms, geometrically similar to the full-scale 
tank (see Fig. 1). The Z-ring stiffeners have a width of ap- 
proximately 3 per cent of the tank diameter and are per- 
forated with lightening holes. The two types of suppressor 
configurations consisted of floating can devices and truncated 
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/ | TANK WALL 


Z-RINGS 


FLOATING CAN 


CONICAL RINGS 
Fig. 2 Suppressor configurations 


Fig. 4 Conical rings mounted in tank 


conical rings, both of perforated material (approximately 25 
per cent open). The cans and conical rings are shown in 
Figs. 2-4. The model suppressor devices were also geometri- 
cally similar to their full-scale counterparts. 


Test Facility 


The excitation system basically involves a rigid platforin 
whose motion can be held to straight line oscillation in a 
horizontal plane, driven by a slider-crank mechanism of 
variable crank-throw and variable crank velocity. Tle 
slider-crank mechanism gives sinusoidal motion to the slider 
element with a maximum error of 2 per cent. The platforin 
oscillates on four pairs of cam-follower rollers, which roll on 
hardened steel bars. Sidewise movement is prevented ly 
four small cam-follower rollers whose axes are at right angles to 
the main rollers. The system is mounted on a heavy test 
frame bolted to a concrete floor, with a dummy platform t» 
provide balance for the system set up at the opposite end cf 
the test frame. The dummy platform oscillates at the same 
amplitude and frequency as the tank platform, but with 1 
phase displacement which produces smooth sinusoidal motion. 
The mass of the dummy platform is matched to that of th: 
tank platform by adding or removing weights. The drive for 
the system consists of a variable speed d-c motor, driven by 
a motor-generator set. 

The complete facility is shown in Fig. 5. It may be men- 
tioned that this facility has been designed and constructe« 
to produce rotational excitation of the model tank, if desired. 


Instrumentation 


Instrumentation is provided to measure the following quan- 
tities: 


1 Tank displacement. 

2 Total force produced by the fluid on tank wall. 

3 Total moment produced by the fluid about tank 
bottom. 

4 Pressure distribution of the fluid on tank wall. 


All quantities are recorded simultaneously on a 16-channel 
oscillograph. The trace of tank displacement, measured by 
the deflection of a thin cantilever beam, is used as a reference 
for phase relationships of all other parameters. 

Total force and moment measurements are considered to- 
gether, since a single system can be designed to provide both 
measurements (7). The measurement is accomplished by use 
of a cantilever beam dynamometer system. The model is 
mounted on four instrumented beams or arms. For force 
measurements, the force strain gages on the individual arms 


Fig. 5 Experimental facility 
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Fig. 6 Comparison of model and ABMA full-scale tests, no damping device, h/d = 0.505 


are oriented in the bridge circuit to be sensitive only to bend- 
ing (or lateral) loads in the arms, and the four arms are 
summed in the bridge to measure the total force. For moment 
measurements, the moment gages on the individual arms are 
oriented in the bridge so as to be sensitive only to axial (or 
vertical) loads in the arms. The outputs from the two forward 
arms in the line of motion are summed in the bridge and out- 
put from the aft arms subtracted. Hence the difference is 
proportional to the pitching moment about an axis in line 
with the pinned platform attachment. 

A major disadvantage to this type of dynamometer system 
is that part of the signal output is due to the mass inertia of 
the tank model, platform and dynamometer arms. These in- 
ertia forces may be of the same order of magnitude as the forces 
due to liquid motion, or in some cases may be even higher; 
thus separate calibration of the inertia response is required. 
The loads resulting from liquid motion are then obtained as 
the difference between two large loads [ (inertia + liquid) — 
inertia]. It is possible, however, to cancel the signals re- 
sulting from inertia loads directly in the electrical circuits. 
The remaining signal is then proportional to the liquid loads 
only, and may be amplified to obtain reliable readings. The 
cancellation is accomplished by providing a separate system 
undergoing the same motion, dynamically similar to the 
dynamometer-tank system but without the liquid loading. 
The output from this “balance’’ system is fed into the dy- 
namometer bridge circuit with a phase difference of 180 deg. 
The balance mass may be adjusted mechanically to nearly 
cancel the inertia signal from the dynamometer system, and 
final trim may be accomplished by electrical attenuation of 
the balance system. A dynamic check of the static calibra- 
tion, using a rigid mass in the tank, established the accuracy 
of the net (remaining) signal response. 

Pressure distributions on the tank wall are measured by 
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means of pressure cells mounted essentially flush with the tank 
wall. Pressure measurements can be made at various heights 
in the tank and along the circumference at one height, for the 
purpose of integrating the pressure distribution for correlation 
with the total force and theory. 

Only the total force measurements shall be discussed in the 
present paper; it is planned to present other measured data in 
subsequent papers. 


Test Results 


Force Data 


Figs. 6-11 compare force data from some of the model tests 
and the ABMA full-scale tests. The liquid height to tank 
diameter ratios of 0.785, 0.595 and 0.505 correspond to the 
height of the liquid above the bottom of an equivalent flat 
bottomed tank of the same volume. In all of these figures, 
the abscissa is the dimensionless excitation frequency dw*/a, 
and the ordinate is the dimensionless total fluid force 
F/(pad*)(X,/d). This form of the dimensionless total force dif- 
fers from that of Equation [9] by the introduction of the 
excitation amplitude, and has been employed only for con- 
venience. Also for convenience, the total force has been given 
in terms of absolute magnitude and phase angle ¢. Theoreti- 
cal response curves applicable, of course, only for the case of 
zero damping, have been shown in some of the figures for com- 
parative purposes. Only a small portion of the test data ac- 
tually obtained is presented in the present paper because of 
space limitations. 

Fig. 6 presents the results of one of a series of tests with no 
damping device.’ The model tests were made with both 


5 The tests referred to here actually do involve the Z-ring 
stiffeners, which may be regarded as a type of damping device. 
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water and methylene bromide. Since the points appear to 
coincide, it is apparent that the Z-ring stiffeners add only 
negligible damping to sloshing characteristics, and therefore 
that viscous effects are negligible for this configuration. This 
is further indicated, in Fig. 6, where the points generally follow 
the theoretical curve corresponding to the case of zero damp- 
ing (8-11).6 Agreement with the ABMA full-scale tests in 
this configuration is very good, except that the full-scale tests 
show considerable scatter. 

Figs. 7, 8, 10 and 11 are comparisons of model and full-scale 
tests with the two types of suppression devices. Since sig- 
nificant damping is present with these configurations and the 
equivalent Reynolds numbers are quite different (1.680 < 10° 
for the model tests and 13.95 X 10° for the full-scale tests), 
the model tests may not be regarded as exact simulation of 
the full-scale tests. Therefore these results should be com- 
pared only on a qualitative basis. Thus, it may be noted 
that although the peak forces in the first mode, comparing 
model and full-scale data, are of the same order of magnitude, 
the large scatter in the full-scale data does not permit any 
accurate comparison of these peak forces. Since the effects of 
damping (except for very large values) are significant only in 
the immediate neighborhood of resonance, it is therefore neces- 
sary to discuss the magnitude of the damping in terms of 
derived damping factors, as given later in the paper. 

Some interesting features of these comparisons, however, 
may be discussed on a qualitative basis. The floating cans 
(Figs. 2 and 3) provided the data of Figs. 7 and 8. Tests 


_® The analyses of these several references are equivalent; the 
theoretical curves shown in the present paper were computed 
from the particular form of the analysis given in (8). 
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were first made with 90 cans (91 were actually used in the 
full-scale tests), which show general agreement with the full- 
scale tests although, again, the latter show considerable 
scatter. Tests were then made with 98 cans, the number re- 
quired to cover the liquid surface completely; these tests 
showed considerably more damping, the effect being more 
pronounced as the liquid height was decreased. An additional 
run was made through the first resonance at h/d = 0.505 
(Fig. 8) with 94 cans; the damping appeared to affect only the 
peak force. 

These tests were made at approximately the same excitation 
amplitude to tank diameter ratio as that of the full-scale tests 
(X,/d = 0.00417 for the model tests vs. 0.00475 for the full- 
scale tests). An additional run at roughly half the previous 
amplitude was made at h/d = 0.505 with 94 cans to investigate 
the effect of excitation amplitude on damping. As indicated 
in Fig. 9, there is little difference in the measured forces except 
at the peak of the first mode. The effect here is considerable. 
The dimensionless peak force at the lower excitation ampli- 
tude is nearly twice that of the higher excitation amplitude, 
showing a large increase in damping with increasing excitation 
amplitude. 

Results of tests with perforated conical rings (Figs. 2 and 
4) are presented in Figs. 10 and 11. As with the 90 floating 
cans, agreement with the full-scale tests is relatively good, 
especially in the region of the first mode, although the full- 
scale test data again show considerable scatter. 

Rather extensive data for tanks in translation or rotation, 
with and without suppressors and applicable to a variety of 
“in-flight” conditions, have been or are currently being ob- 
tained as part of the present program, and may be presented 
in subsequent papers. 
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Fig. 12 Variation of damping factor with equivalent Reynolds number, Xo/d = 0.0045 


Damping Factors 


The damping provided by these two types of suppression 
devices was measured from the force response data, from 
tests with two sizes of models (14.4 and 25.0 in. in diameter) 
by an iterative procedure based on the half bandwidth tech- 
nique applied to the test data, as given and described in (12). 

Fig. 12 shows a plot of damping factor g vs. the effective 
Reynolds number, using values of density and viscosity at 
room temperature taken from published tables (6).’ It is 
noted from this figure that simulation of full-scale conditions 
is easily possible for kerosene at 1 or more G by employing a 
slightly larger model tank, but that simulation of full-scale 
tests for water at 1 G is not possible with any of the model 
fluids employed. Values of pd’/*a’/*/u for the various model 
tank-liquid combinations are given in Table 1. 

It may be further noted from Fig. 12 that the damping due 
to the 90 floating cans is consistently lower than that due to 
the conical rings. The damping factors at pd’/a'/?/u = 
13.95 X 108, corresponding to the full-scale water tests made 
at ABMA, and also to the same tank with kerosene at an ac- 
celeration of 10 G, are obtained by extrapolation. These data 
on damping factors appear to be relatively consistent, and it is 
believed that the extrapolation is reasonable. Comparisons 
of damping factors, as obtained by this extrapolation from 
model data (Fig. 12), and as obtained from iteration of force 
response and comparison with a family of theoretical curves 
of the ABMA tests (12), are given in Table 2. 

Damping‘factors obtained from the model tests with both 
floating cans and conical rings agree reasonably well with 
those obtained from a comparison between the full-scale tests 
and a family of theoretical curves. The damping factors ob- 
tained direetly by iteration from the full-scale tests are not 
considered to be reliable because of the large scatter of the 
data, and are given here only to complete the comparison. 


’ The viscosity of methylene bromide given in (6) is incorrect 
and should be multiplied by 10. This error was pointed out to the 
authors by Professor Paul E. Sandorff of MIT. 
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The former procedure introduces some smoothing of the data 
and should provide more accurate values, although there is 
some evidence to indicate that the damping factors thereby 
obtained are somewhat high. Thus, a reduction of those 
values would tend to improve even further the agreement with 
those obtained from the model test data. 

The essential result to be noted from Fig. 12 is, however, 
that the damping factors are quite definitely dependent upon 
the equivalent Reynolds number. Although there is some 
scatter in the damping factors obtained from the model test 
data, this does not negate the principal feature involved. 
These results demonstrate quite clearly that model studies in 
which damping or slosh suppression devices are involved can- 
not be conducted without due regard for the requirements 
imposed by similitude theory. 
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Effective Reynolds numbers for various models 


Effective Reynolds number, 
pd?/2q'/2/p 


Table 1 


Model tank—liquid 


small—water 0.686 < 10° 
large—water 1.565 10° 
small—methylene bromide 1.680 x 10° 
small—methylene chloride 1.970 x 108 

4.500 10° 


chloride 


Table 2 Damping factcr gy at pd*/2a'/?/u = 13.95 X 106 | 
——ABMA Results 


SwRI results 


From compari- 
son with theo- 


By iteration 
from full-scale 


from Fig. 12 _ retical curves experiment 
h/d cans rings cans rings cans rings 
0.505 0.125 0.200 0.136 0.28 0.20 0.25 
0.595 0.125 0.175 0.146 0.28 0.25 0.28 
0.785: 0.125. 0.152 6:16 0.21 
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Conclusions 


The data presented in the present paper permit several 
conclusions to be drawn: 

1 Model experiments performed in accordance with the 
similitude requirements derived herein can be used to provide 
quantitative data which can be applied directly to full-scale 
tanks, even with large amounts of damping present. The 
relatively good agreement obtained between the model and 
full-scale tests lends confidence to the prospect of simulation 
of “in-flight”’ sloshing characteristics. 

2 Stiffener rings with lightening holes, of the order of 
width used in this configuration (3 per cent of tank diameter), 
add negligible damping to sloshing forces. 

3 The amount of damping produced by both floating cans 
and conical rings is strongly dependent upon excitation 
amplitude. For this configuration, the perforated conical 
rings provide more damping than the 90 floating cans, and 
approximately the same amount of damping as 94 cans (at an 
excitation amplitude of Xo/d = 0.00417). The damping pro- 
vided by conical rings is furthermore dependent on the depth 
of liquid in the tank. 
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Rocket Boost Trajectories for 


Maximum Burnout Velocity 
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The problem considered is the definition of optimal trajectory programs which allow a given 
rocket boost vehicle to launch its payload with maximum velocity, for prescribed burnout values 
of the path angle and altitude. Methods of the calculus of variations are employed to achieve the 
desired solutions. The equations of motion consider a spherical nonrotating Earth model and in- 
clude the effects of aerodynamic forces. Inequality constraints suitable for satisfying certain 
practical vehicle design and physiological limits are introduced. An optimal coasting period of 
initially unspecified duration may be included in the trajectory to provide increased burnout alti- 
tude. Numerical results illustrating the utility and validity of the method are included. 


HE ADVENT of spaceflight and quasi-spaceflight places 
particularly stringent requirements on the performance of 
rocket boost systems. Maximum mission performance is 
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realized only by insuring the most efficient transfer of impulse 
from the given boost system to the spacecraft. This goal 
may be achieved by properly programming the boost tra- 
jectory, with due consideration of the constraints imposed 
by the vehicle design and the capabilities of the control system 
employed. It is currently common practice to program boost 
trajectories by what may best be described as “educated” 
guessing of the values of certain control parameters. 

This paper presents a method for the determination of 
optimal boost trajectory programs (extremals) which yiel! 
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the maximum burnout velocity for a specified booster, pay- 
load and mission. The proposed method is useful in the pre- 
liminary booster design period for defining ultimate per- 
formance capability, indicating efficient means of utilizing 
the available control systems and providing a systematic 
basis for evaluating performance trades resulting from design 
modifications. 

Specifically, the proposed method yields the time history 
of the two-dimensional boost trajectory which results in the 
maximum burnout velocity attainable by a given multistage 
booster configuration, subject to specified initial and final 
boundary conditions. Initial conditions used (velocity, path 
ang! and altitude) allow for either ground or air launching; 
the final conditions include the path angle and altitude ap- 
propriate to the mission under consideration. Effects of 
aero\lynamic forces are included. The trajectories may be 
constrained to satisfy limits imposed by structural deforma- 
tion. longitudinal stability, aerodynamic heating and acceler- 
ations acting on the crew of the spacecraft. An optimal 
coasting period may be included in the extremal, if desired. 

Solutions derived herein are based upon application of 
metiiods from the calculus of variations. As such, the gen- 
eral approach to the problem is similar to those used by Lieb- 
hol’ and Brown (1) and Leitmann (2), who were concerned 
wit): the derivation of boost trajectories yielding maximum 
rane for ballistic missiles. Inclusion of an optimal coasting 
stage between powered stages is made possible by a special 
interpretation of the basic variational problem. 

Use of a high speed digital computer is required to calculate 
numerical solutions by the proposed method. In the interest 
of promoting efficient computer utilization, a rapidly con- 
verging iteration technique is adapted to the proposed 
method. 

Two numerical examples are included to illustrate the 
validity and utility of the proposed method. In these ex- 
amples, extremals are determined for the launching of circular 
satellite and boost-glide vehicles by a hypothetical booster 
representative of a moderately advanced solid rocket. 


Equations of Motion 


In this analysis the rocket is taken to be a point mass mov- 
ing in the vicinity of a spherical nonrotating Earth (or other 
celestial body). Orientation of the rocket relative to Earth 
and the force, velocity and acceleration vectors which exist 
during flight are illustrated in Fig. 1. The trajectory of the 
rocket satisfies the following equations 


md = T cosa — D — mg sin y 


md - = Tsina + L — mgcos y 


where the dots above the symbols denote time derivatives. 
The first two of these are the equations of motion for tra- 
jectory oriented reference axes whereas the third is a purely 
kinematic relationship required to allow the use of four de- 
pendent variables (#, y, h and qa) in the problem. Since 
the time derivative for the angle of attack does not appear in 
Equations [1], the angle of attack @ is the control variable 
(i.e., driving function) in the problem. 

A spherical nonrotating Earth model was selected as a com- 
promise which provides useful trajectory results without 
introducing the complications of oblateness and rotation. 
Approximate gross effects of Earth’s rotation may be applied 
to solutions derived by the proposed method for closer esti- 
mates of burnout conditions. Effects of thrust deflection 


* Numbers in parentheses indicate References at end of paper. 
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ARROWS DENOTE 
POSITIVE DIRECTION 


Fig. 1 Orientation and vector system of rocket 


were also ignored, principally to eliminate an additional 
dependent variable from the problem. In most instances the 
errors resulting from neglecting thrust deflection will be well 
within the accuracy levels ordinarily associated with pre- 
liminary design studies. 

For the purposes of this analysis it is assumed that the 
mass, thrust, lift, drag and acceleration of gravity are known 
at all points of the boost trajectory and have the following 
prescribed functional dependence: m = m(t), T = T(h, 0), 
L = L (8, a, h), D = D(d, a, h) andg=g(R+h). Note 
that this excludes the effects of aerodynamic lag in the appli- 
cation of the aerodynamic forces. 


Trajectory Optimization 


Prior to initiating the analytical development of an opti- 
mization technique, due consideration must be given to the 
selection of the parameter to be extremized. This involves 
the character of the vehicle considered as well as its general 
mission capabilities. Numerous papers appearing in the 
literature discuss the single mission concept, i.e., maximizing 
such things as ballistic range, satellite orbit altitude, burnout 
height, etc., leading the uninitiated reader to believe that 
variational techniques are somewhat specialized. The value 
of any vehicle/mission optimization technique is directly 
dependent on the overall utility of the technique in designing 
vehicles or programming their trajectories for a variety of 
missions. A primary goal of the subject study involves 
formulation of a technique having wide utility. 

The principal purpose of a rocket booster is to convert the 
latent energy of its propellants into useful energy of the pay- 
load. In this sense, useful energy is defined as the energy 
required to perform a given mission most efficiently. It is 
not to be confused with maximum energy which is, in most 
instances, neither required nor desired. For sensible mission 
planning, the useful energy so imparted must be suitably 
directed by specifying the direction of the velocity vector and 
the burnout altitude. In this format, the proposed tech- 
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nique to be described provides valuable design information 
for the many types of missions whose extremals are members, 
defined by appropriately specified values of burnout altitude 
and path angle, of the family of maximum velocity tra- 
jectories. 

Rocket burnout velocity is determined by integration of the 
first of Equations [1] 


= +o [2] 


The problem considered here is to determine those functions 
3(t), y(t), A(t) and a(t) which maximize the integral [2] sub- 
ject to certain specified boundary conditions and constraints. 
Initial boundary conditions used specify the velocity, path 
angle and altitude at the beginning of the trajectory, ic. 


at t = b, h=h [3] 


Terminal boundary conditions used are: 
att = ty, h=H [4] 


Values of T and H are deduced from the mission being con- 
sidered. 

Three of the constraints employed require that Equations 
[1] be satisfied at all points of the trajectory. These equa- 
tions are written in the form 


= 0 — FY, y,h,a,t) = 
= — GY, y, h,a,t) = 
=h—Usiny =0 


oo 


Solutions obtained with use of only these three constraints 
are of interest, but may be impractical when the physical 
limitations of the vehicle or its crew are considered. Such a 
contingency is anticipated, and may be avoided, by the 
addition of an inequality constraint of the form 


Gs = Ln — In(8, Y, h, > 0 m=1,2,3... [6] 


where each /» must be an explicit function of a, and L,, is the 
corresponding limiting value. 

Since there are four dependent variables in the problem, 
using the equality form of just one of the constraints [6] 
together with Equations [5] uniquely determines the boost 
trajectory. For practical applications a number of limita- 
tions are prescribed for a given problem. In the event that 
more than one limit is violated at a given time, the com- 
puter selects and uses the most critical. An extremal must 
be free of all the constraints [6] for some portion of the 
powered flight. Otherwise, the variational problem degener- 
ates into an ordinary assumed trajectory program problem. 


Design and Physiological Constraints 


Before proceeding with further discussion of the variational 
problem it is appropriate to note some types of limitations 
which may be satisfied by use of the inequality constraints 
[6]. One of the more obvious vehicle design limitations is 
the normal load factor. This may be satisfied by a dual 
purpose constraint. Let the angle B be measured positive 
aft of the normal to the vehicle longitudinal axis. Then, 
the acceleration vector acting in the 6-direction has the mag- 
nitude 


m m m 


When 8 equals zero deg. Equation [7] gives the normal acceler- 
ation. 

In manned vehicles, the principal physiological limitations 
during boost are prescribed by human tolerance to accelera- 
tion, which is affected significantly by the direction of the 
acceleration vector. Thus, Equation [7] may be used to 


614 


satisfy such limitations, provided that 8 is specified in the 
direction of the least human tolerance to acceleration. 

Other vehicle design limitations are imposed by structural 
deformation, longitudinal stability and aerodynamic heating. 
Accurate definition of these limits requires complex analyses 
which are not appropriate to the present study. However, 
trajectory constraint may be achieved without undue com- 
plication by use of certain similarity-type parameters. 

A structural bending of the boost vehicle may be caused by 
aerodynamic loads which, for most boost trajectories, are 
critical at relatively low altitudes and moderate supersonic 
speeds. If the booster is aerodynamically stabilized or is 
carrying a winged payload, the aerodynamic loading acting 
on the fins or wing produces most of the bending moment. 
This loading is nearly proportional to the product of the 
dynamic pressure and angle of attack. Hence, the relation- 
ship 


= (1/2)pd? | [8] 


may be used to control the structural bending. 

Vehicle static longitudinal stability is ordinarily acceptable 
if the angle of attack remains within certain limits at a given 
Mach number. Thus, an appropriate control parameter is 


= |a| 


In this case L3 may be specified as a function of Mach number. 
Aerodynamic heating limitations may be satisfied by using 
a control parameter of the form 


ly = Tev(8, h, a) [10] 


where 7. is the equilibrium wall temperature at an adroit!y 
selected point of the vehicle’s surface. The presumption here 
is that the equilibrium wall temperature at the selected 
point will be representative of that over the remainder of 
the body and, therefore, serves as a suitable control over the 
total heat absorbed. Forms of equations used for computing 
T will depend on the vehicle surface shape and the location 
selected for the monitored point. A general representation 
may be achieved by approximating 7’. in Equation [10] 
with a polynomial surface within a critical region of the #, 
h, aspace. 

The control functions noted in Equations [7 through 10] 
are used solely to modify unsuitable extremals. That is, an 
unconstrained extremal is determined first and the controls 
are used only as required to derive a constrained optimal path 
which satisfies the pertinent vehicle design and physiological 
limitations. 


Euler-Lagrange Equations 


The variational problem is stated as a modification of 
Lagrange’s form (3); i.e., it is required to find a stationary 
value of the integral 


[= (9 4 es) dt f, (0 + her) dr 


where the Lagrangian multipliers ; are time dependent and 
i, represents the variable time at the end of the optimal coast- 
ing period. Note that the relationship Aug; = 0 must be 
satisfied in this formulation. A necessary condition (proof 
of sufficiency is beyond the scope of this paper) for the exist- 
ence of such a stationary value is that the first variation of 
J must vanish 


+ f, 6F*(r) dr =0 [12] 
where 


F*(t) =x + Xk, 
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with the notation xx, k = 1, 2, 3, 4 representing 3, y, h and 
a, respectively, being adopted for brevity. In the first in- 
tegral of Equation [12] the trajectory comparisons are made 


at the same values of time. Hence 
4 oF* 
ka 
Furthermore 


h OF* 
— — (6xx — xx6t)dt = 


oF* 
(6x — xan | - (6x — — at [14] 
UXk 


In the second integral of Equation [12] the comparisons are 
made at the same values of 7, and results corresponding to 
Equations [13 and 14] may be derived by replacing ¢ with 7 
in those equations. Substituting these results in Equation 
{12| yields the complete expression for the first variation 


oF* 


OXe 

oF 


Since the variations 6x; are arbitrary, application of the basic 
lemma of the calculus of variations yields 


oF* doF* 


dt 


k = 1,2,3,4 16 
Ox~ dt OX: [16] 
or, in terms of the physical variables 
oF oG 
Moy 
or oG 
— 2— =0 
"da * da oa 


These are the Euler-Lagrange equations for the problem, 
which make both integrals of Equation [15] vanish. 

Inspection of Equations [17] reveals that they are homo- 
geneous in the multipliers \;. This property, which was 
also noted by Cicala and Miele (6) and Breakwell (4), makes 
solutions of these equations independent of the initial values 
of one of the multipliers. In the present analysis, A; was the 
multiplier selected. The Euler-Lagrange equations retain the 
same form except that the multipliers and their time deriva- 
tives are replaced by 


= di [18] 


The initial values A), are arbitrary. 


Corner Conditions 


Staging points of the trajectory are characterized by dis- 
continuities of arbitrary magnitude in the thrust and mass 
of the rocket. In the physical case (and in the variational 
problem) it is required that 3, y, h and h be continuous at 
staging, whereas 3 and y may be discontinuous. Conclu- 
sions regarding some of the other parameters are deduced 
from the Erdmann-Weierstrass corner conditions. 


1960 


* 
— x,0t)dt +f, (= - 


Integrating the Euler-Lagrange equations, [16] yields the 
expressions 


ax k 


+C 


[19] 


k = 1,2,3 


where C is an integration constant. Since F* is bounded 
and continuous within each stage of the trajectory, it follows 
that the values of \;(¢ = 1, 2, 3) are continuous over the entire 
trajectory, including staging points. 

The corner conditions yield no information about a@ and 
As. When the trajectory is not constrained by one of the 
Lm limits (¢, > 0 and A; = 0), the fourth Euler-Lagrange 
equation of [17] shows that @ may be discontinuous at 
staging. This is the result of assuming the rocket to be a 
point mass without inertia and including aerodynamic forces 
in the analysis. Though an instantaneous jump in angle of F 
attack is unrealistic, its effect on the overall trajectory is in- 
significant. When the trajectory is constrained by one of 


4 ort 


OX: 


d 
dt aX: 


OX: 


the L» limits (g, = Oand A, ¥ 0), amay be either continuous 
or discontinuous at staging, depending on the form of control 
parameter J, which is involved. The multiplier 4; may also 
be discontinuous at staging. In view of these circumstances, 
the character of the trajectory may change at staging, i.e., 
from constrained to unconstrained. 


End Conditions 


Satisfying the remainder of Equation [15] provides the re- 
quired end conditions for the problem. Applying the corner 
conditions and noting that no initial variations are allowed, 
Equation [15] becomes 


G Xk n+ ox) 


where ¢. is the total coasting time and 67 = 0, since the burn- 
ing times of the stages following the coasting period are speci- 
fied. By definition, \;¢,; equals zero on trajectories, so Kqua- 
tion [20] may be written as 


=0 [20] 


—(Md + Avy + Ste + + + = 0 
[21] 


Final values of y and h are specified by conditions [4], so 
6y, = 6h; = 0. Furthermore, on an extremal, 6%; is equal 
to zero, so Equation [21] reduces to 


(Aig ey + Ash) = 0 [22] 


for extremals containing an optimal coasting period. Hence, 
Equations [4 and 22] define the desired extremal in this case. ' 
When such a coasting period is not included, 6¢, of Equation 

[21] is also zero and the only conditions required are those 

given in [4]. In the latter case every trajectory computed 

is an extremal for the particular end conditions achieved. 


Extremal Solutions 


In the preceding section the problem under consideration 
has been stated in terms of the initial conditions [3], the ter- 
minal conditions [4] (and [22], if desired) and the set of Equa- 
tions [5,6 and 17]. Thus, it is seen that the boundary value 
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problem involved has its boundary conditions applied both 
at the beginning and end of the trajectory. - The Lagrangian 
multipliers provide, in effect, dummy initial conditions which 
permit the initiation of a given trajectory integration. 
Consider the problem of finding an extremal which con- 
tains a coasting stage. In this case the appropriate equations 
of the set [5, 6 and 17] are solved, using the initial conditions 
[3] in conjunction with guessed values of A»*, Aa* and ¢., to 
determine a given trajectory. The desired extremal is found 
by iterating on the guessed values until the terminal condi- 
tions [4 and 22] are satisfied within prescribed tolerances. 
When a coasting stage is not included, the problem is handled 
in an analogous fashion, with ¢, = 0 and condition [22] 
eliminated. Methods used for trajectory integrations and 
iteration for the desired extremal are next described. 


Trajectory Integrations 


All trajectory integrations are performed by means of the 
Runge-Kutta four point scheme for several dependent vari- 
ables, using time as the independent variable and beginning 
with the given initial conditions. Trajectory segments which 
are computed fall into four categories: Unconstrained, con- 
strained, mixed and initial. 

Unconstrained trajectories are true variational paths which 
are not constrained by any of the L» limits. In this case 
X, is equal to zero and the integration involves the simul- 
taneous solution of equations, [5 and 17]. The angle of 
attack at any instant is calculated by an iterative solution 
of the last equation of [17]. 

Constrained trajectories are defined as those constrained 
by one of the prescribed Lm limits. The path is uniquely 
determined by the solution cf Equations [5 and 6], with a 
being determined from Equation [6], either directly or by 
iteration. The Euler-Lagrange equations, [17] are then 
solved independently, with \, ~ 0, to determine the values 
of the Lagrangian multipliers. 

Mixed trajectories have both unconstrained and con- 
strained segments. Computations for these segments use 
the techniques noted, but special tests are required to ac- 
commodate entering and leaving a constrained segment. 
Locating the point where a trajectory first becomes con- 
strained involves an interpolation to determine the time 
when Equation [6] becomes an equality. A similar inter- 
polation is employed to determine the time when the tra- 
jectory leaves a constrained segment, except that A, = 0 
must be satisfied in this instance. When more than one 
limit is imposed in a given problem these tests are coupled in 
a manner appropriate for selecting the most critical. 

When the velocity is zero the set of equations involved in 
the proposed method becomes singular. Hence, for ground 
launched rockets, an initial, pre-programmed trajectory 
segment spanning the time interval 0 < ¢ < f& must be com 
puted by solving either y = Yo, or y + @ = constant, together 
with Equations [5]. Lagrangian multipliers are not com- 
puted during this period, and the variational problem be- 
gins at t = &, with v = v(t), Yo = Y(t) and ho = h(t). The 
angle of attack will ordinarily be discontinuous at this time. 


Iteration Procedure 


Use of an efficient iteration technique is mandatory for 
finding desired extremals with reasonable expenditures of 
computer time. The technique adopted for use with the pro- 
posed method was attributed by Montague (5) to Mark 
Robinson. To illustrate the adaptation of this technique, 
write the desired end conditions for an extremal with a coast- 


ing stage as 


yn=H-h; =0 


f= (ud + dey + Ash) mt, =0 


The procedure begins with the computation of four trajecto- 
ries having the following initial conditions and final results 


+ Arz,*, te > My 
aa, ag + Ad;,*, 
The first of these trajectories defines the residuals produced 
by the guessed values of d2,*, As,* and ¢&: R, = &, 


R, = », and R; = ¢,, and the subsequent three trajectories 
define the derivatives 


og & & on — 


§ 


etc. To find the desired extremal, it is necessary that al of 
the residuals vanish. For assumed linear variations betw: en 
the initial and final conditions, this will be true if, in matrix 
notation 


|} = r=1,2,3 s=1,2,3 [24] 


where pi, p2 and p: are A»,*, \3,* and ¢., respectively, and 
on 
Oy* dle 
ar ak 


Solving Equations [24] yields the next approximations to the 
initial values 


|Acr| 


Pr’ = pr + Ap,’ 
that are used to compute the fifth trajectory. 

This procedure, which may be cycled to accomplish the 
complete iteration, has the disadvantage of requiring three 
redundant trajectory computations per correction. Such 
redundancy is avoided by using the matrix modification 
procedure suggested by Robinson. The modified matrix 
satisfies the relationship 


+ = frre) — Raat 


where the AA,, represent modifications of the individual 
matrix elements, and the residuals R,(p,’) are computed from 
results of the fifth trajectory. Subtracting Equation [25] from 
[24] yields 


= 
which has a solution (non-unique) of the form 
AAs = [R.(pr’)| Ap,’ / (Ap,’)? 
The next corrections are 
Ap,” = |Ag + 
and the next approximations to the initial values are 


pr” = Pr’ Ap," 


Subsequent approximations are computed in an analogous 
manner, with one trajectory integration per correction. The 
same iteration technique, with appropriate reduction to two 
parameters, is used for finding extremals that do not have 
coasting periods. 


Numerical Results 


The method previously described was coded on an IBM 
704 computer at Convair-Fort Worth and is now being use 
in design studies. To confirm the analytical formulation 
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Fig. 2 Extremum selection for boost-glide vehicle launching, 
r = 0, Wy = 4000 Ib 


employed, numerical results for two hypothetical missions 
are presented and analyzed below. In these calculations the 
mass flow rate and thrust were represented by 


m = —(1/e)(TR + AezPr) [26] 
and F 
T = Tr + Aer(Pr — P) [27] 
where 
Tp = rated thrust (may be time variant) at an altitude hav- 
ing the ambient pressure Pz 
Apr = nozzle exit area 
c = effective exhaust velocity 
P = ambient pressure at the operating altitude 


The general representations used for the lift and drag coeffi- 
cients are 
Cy = (Ci, + K| sin @|) sina [28] 


Cp = Cp, + Cx sin a [29] 


In the examples next discussed, K = 0. The particular forms 
[28 and 29] were selected to make the fourth Euler-Lagrange 
equation of [17] be continuous for all @ values. 


Boost Rocket 


Two missions were selected for the illustrative examples. 
These consist of launching a hypersonic glider and a satellite, 
each weighing 4000 Ib from the surface of Earth. Charac- 
teristics of the boost rockets employed are given in Table 1. 


Table 1 Boost rockets 
Stage 1 Stage 2 Stage 3 
rated thrust, 
lb X 1073 840, sea lev. 140, vacuum 15.12, vacuum 
specific impulse, 
sec 250, sea lev. 280, vacuum 280, vacuum 
burning time, sec 60 60 60 
stage weight, 
lb X 10-3 219.1 32.6 3.6 
propellant weight, 
Ib X 1073 201.6 30.0 3.24 
nozzle exit area, ft? 25 4 re 
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These hypothetical vehicles are representative of current 
solid rocket technology. Stages 1 and 2 were sized to ac- 
commodate the boost-glide mission, and stage 3 was subse- 
quently added to achieve a low satellite orbit. In the ex- 
amples, the rated thrust was assumed constant and realistic 
values of the lift and drag coefficients were employed in con- 
junction with the ARDC 1956 atmosphere model. 


Trajectory Selection 


Selection of the desired boost trajectories simply consists 
of matching the booster capability with the general mission 
requirements. Fig. 2 illustrates the mode of trajectory selec- 
tion employed for the launching of a hypersonic glider. The 
broken lines in this figure represent equilibrium glide paths 
computed by specifying y = 0 in the second equation of [1]. 
The two-stage booster capability, labeled the locus of ex- 
tremals, was computed by specifying I = 0 and several 
burnout altitudes, with the initial pitch point arbitrarily set 
at 10 sec. Some difficulty was encountered in ascertaining 
suitable values of \»,* and A;,* to achieve convergence to the 
desired extremal in these calculations. This difficulty was 
attributed to the sensitivity of the trajectories to the initial 
values of the Lagrangian multipliers, caused by specifying 
burnout altitudes in regions where the aerodynamic forces 
are significant. 

One of the computed trajectories has the time history indi- 
cated in Fig. 3. End conditions of this trajectory are indi- 
cated by the circled point of Fig. 2, having h = 250,000 ft, 
v = 24,653 fps and W/SC, = 288. The relatively short 
burning time of the rocket caused considerable turning 
throughout the entire trajectory. 

Trajectory selection for the satellite launch example is 
illustrated in Fig. 4. In this instance the required condi- 
tions are given by the circular satellite velocity relationship 


= + hy) 


where the values of R and the sea-level acceleration of gravity 
were those of the ARDC 1956 atmosphere model. The 
locus of extremals was computed for T = 0 and several values 
of H, with the coasting time optimization included. In this 
case no convergence difficulties were encountered, with the 
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desired extremal being found after an average of 12 trajectory 
integrations, a small number for a three-parameter iteration. 
The intersection of the circular satellite velocity line with 
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Fig.4 Extremum selection for satellite launching, ' = 0, W, = 
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the locus of extremals represents the maximum orbit altitude 
attainable (227 nautical miles) with the particular rocket- 


payload configuration assumed. Fig. 5 shows the time 
history of the maximum orbit altitude trajectory, having 
an optimum coasting time of 416 sec. This figure indicates 
that the gravity turn program, highly touted for liquid 
rockets, is hardly optimum for solid rockets of the type as- 
sumed in this example. The short burning time of the first 
two stages does not permit sufficient gravity turning prior to 
coast. 

Both of the trajectories indicated in Figs. 3 and 5 could be 
closely simulated by realistic guidance and control prograins. 
For example, three constant attitude angle rates (y + 4), 
one in each powered stage, would yield very nearly the same 
burnout conditions. Losses in final velocity caused by :p- 
proximating the optimum control program will be small siice 
the optimization technique requires that the derivative of 
final velocity with respect to instantaneous angle of atteck 
be equal to zero. 


Summary of Results 


An evaluation of the numerical results shown is summarized 
in Table 2, where the ideal velocity—velocity losses caused 
by gravity, turning and drag—and fraction of ideal velocity 
attained are shown, respectively, from left to right. Ideal 
velocities shown include the effects of thrust variations caused 
by ambient pressure changes along the respective trajectory. 
The high percentages of ideal velocity obtained substantiate 
the validity and utility of the method. 


Table 2 Evaluation of numerical results 
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Fig. 5 Time history of extremal for satellite launching 


Per 
cent | 
Bid Ad, Ava Adp Via 
Glider launching | 
stage 1 13,085 —1,458 -—125 —-1,345 77.6 
stage2 15,430 —314 -415 -—205 94.0 | 
mission 28,515 -—1,772 -—540 —1,550 86.5 
Satellite launching 
stage 1 12,655 —1,464 -—1,275 77.7 
stage 2 12,350 —527 —30 —80 94.8 | 
coast — 1,463 
stage 3 5,005 —1 100.0 | 
mission 30,010 -—3,455 —-110 —1,355 83.6% | 


* Including coast stage loss. 
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Oblique Shock Relations for 
Air at Mach 7.8 and 7200 R 
Stagnation Temperature 


H. T. NAGAMATSU,? 
J. B. WORKMAN*® and 
R. E. SHEER Jr.* 


Research Laboratory, 
General Electric Co. 
Schenectady, N. Y. 


Oblique shock relations for dissociated air in thermodynamic equilibrium have been calculated 
for a free stream flow Mach number of 7.8 in the test section with an equilibrium stagnation tem- 
perature of 7200 R in the reservoir. The results for the density, pressure and temperature ratios, 
and flow deflection across the shock wave are presented as functions of the shock wave angle. Com- 
plete equilibrium has been assumed in the calculations, utilizing the best available thermodynamic 
properties of air for the region considered. The real gas oblique shock relations have been ex- 
perimentally verified by testing an adjustable two-dimensional wedge model in a hypersonic shock 
tunnel. Correlations of the calculated and experimental pressure ratios and shock wave angles are 
presented as functions of flow deflection angle. When allowance was made for the boundary layer 
displacement effect, the correlation was seen to be good. The reservoir conditions were selected to 
insure that the flow in the shock tunnel would be close to equilibrium. 


fen ADVENT of missiles and space vehicles traveling at 
hypersonic flight speeds has recently stimulated much in- 
terest in the shock waves produced by such bodies. In par- 
ticular, there has been considerable discussion of the “real gas 
effects” in the air immediately behind high Mach number 
shock waves. By “real gas effects’ is meant the departure of 
the thermodynamic properties of the medium from a perfect 
gas with a constant specific heat ratio of 1.4. Although air 
departs from this perfect gas model at temperatures as low 
as 540 R, this concept has been used with considerable 
success in the past for analyzing low Mach number shock prop- 
erties. However, at any flow condition even approaching 
that under consideration here, this model must be abandoned 
and, instead, the actual thermodynamic properties of air used. 

As soon as shock waves strong enough to excite vibration, 
dissociation or ionization and strong enough to produce chemi- 
cal reactions are considered, the reaction rate question im- 
mediately comes to the forefront. In dealing with shock 
processes wherein a gas is required to absorb large quantities 
of kinetic energy as internal energy virtually instantaneously, 
one must be concerned not only with the final thermodynamic 
properties, but also with the time needed to achieve thermo- 
dynamic equilibrium. A rigorous analytical method for 
treating the shock problem would be to set up a careful 
kinetic theory model of the gas molecules with reaction rates, 
and then to analyze exactly what happens as the flow proceeds 
through the momentum and energy readjustments at the 
shock region. Unfortunately, there is insufficient knowledge 
of the properties and reaction rates of air at the present time 
to employ such a method. It is conceivable also that even if 
all of the rate constants and interaction parameters were 
known, the complexity of the calculation would prevent a 
practical solution. 
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Recently (1,2)5 shock relations for air under the assumption 
of instantaneous thermodynamic equilibrium downstream of 
the shock wave for two-dimensional conditions have been 
calculated. But, only a negligible amount of experimental 
data is available to correlate with these calculations. This 
paper presents the oblique shock relations for the assumption 
of thermodynamic equilibrium downstream of the shock 
wave, using the high temperature thermodynamic properties 
of air given in (3 and 4). The calculated shock wave angles 
and surface pressures are compared with the experimental 
data obtained in a hypersonic shock tunnel (5) on an adjust- 
able wedge model. 


Analysis of Oblique Shock Relations in Real Gas 


The thermodynamic properties of air across the range of 
shock wave angles have been calculated using a procedure 
similar to that discussed by Moeckel (1) and others. The 
basic concepts of conservation of mass, momentum and energy 
were applied assuming the fluid to be a continuum, inviscid 
and always in thermodynamic equilibrium. If it is assumed 
that equilibrium is attained immediately behind the shock 
wave, the wave may be assumed to be infinitely thin for pur- 
poses of analysis. In solving for the conditions after the 
oblique shock wave, the available references on equilibrium 
thermodynamic properties for air were utilized. For tempera- 
tures above 3600 R the information presented by Feldman (3) 
was employed; from 540 to 3600 R, Hirschfelder and Curtiss 
(4). Below 540 R the air was assumed to be a perfect gas with 
a specific heat ratio of 1.4. 

The oblique shock wave relations for air in thermodynamic 
equilibrium were calculated for the following free stream 
condition: A flow Mach number of 7.8, ambient temperature 
of 810 R, static pressure of 0.0206 psia, and stagnation tem- 
perature after a normal shock wave of 6070 R. These values 
were obtained in the test section of a hypersonic shock tunnel 
(6) with a reservoir equilibrium stagnation temperature of 
7200 R, a pressure of 500 psi and a nozzle area ratio of 576. 


5 Numbers in parentheses indicate References at end of paper. 
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For this particular test condition, both static and impact 
pressures were measured along the axis of the conical nozzle. 
The results agreed very closely with the calculated values ob- 
tained assuming equilibrium during the expansion process. 
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Fig. 1 Calculated Mach number behind shock and flow deflec- 
tion angle vs. shock wave angle 
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Fig. 2 Calculated temperature ratio across shock vs. shock 

wave angle 
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Fig. 3 Calculated pressure ratio and density ratio across shock 
vs. shock wave angle 


The relation between the flow deflection angle across the 
shock wave @ and the shock wave angle @ is shown in Fig. | 
for these initial conditions. For shock wave angles greater 
than 30 deg the flow deflection for y = 1.4 is much less than 
that for the real gas. The maximum flow deflection angle at 
detachment is 54.5 deg for the present case, whereas for a con- 
stant specific heat ratio of 1.4 and a flow Mach number of 7.8 
the deflection angle is only 43.5 deg. This large difference in 
the flow deflection angle is caused by the increase in the 
density ratio across the shock wave and by the decrease in 
the ratio of specific heats which occurs at these elevated 
temperatures. 

For the same upstream flow Mach number of 7.8, the flow 
Mach number after the shock wave is always less for the 7 = 
1.4 case than for the real gas case, up to a nearly normal shc ck 
wave as indicated in Fig. 1. It is also apparent from t iis 
figure that for a Mach number of unity the shock wave ange, 
after the shock, has increased from 65 to 72.5 deg. Thus, with 
a detached shock wave for a symmetrical blunt body at hyp: r- 
sonic flight conditions, the sonic point on the shock wave \ ill 
move toward the axis of symmetry as the flow departs frm 
y = 1.4. The flow Mach number after a normal shock wa\e, 
B = 90 deg, decreases from 0.39 to 0.30 owing to the rial 


gas effects. 
In Fig. 2 the temperature downstream of the shock is 
shown as a function of the shock wave angle for both y = !.4 


and the real gas conditions. For the same flow Mach nuin- 
ber, the perfect gas temperatures are higher than those for tie 
real gas flow. Because of the energy absorbed in dissociation, 
chemical formation and ionization, the equilibrium tempera- 
ture is lower. The greatest temperature difference between 
the two cases occurs when the shock wave is normal to the 
free stream flow. For an infinitesimal disturbance in the free 
stream condition, with the shock wave angle equal to the 
Mach angle, the temperature ratio for both cases is the same, 
since there is no real gas effect. But as the shock wave angle 
becomes greater, corresponding to a stronger shock, the tem- 
perature after the shock increases and the real gas effect. be- 
comes greater. Thus, with an adjustable wedge it is possible 
to have the free stream condition ahead of the shock wave to 
be that of a perfect gas, and after the shock wave to control 
the degree of real gas effect by adjusting the flow deflection of 
the wedge angle. 

The pressure and density ratios are shown in Fig. 3 as func- 
tions of the shock wave angle for both the perfect gas, y = 
1.4, and the real gas case. The real gas effects increased the 
pressure after the strong shock waves. This corresponds to 
large shock wave angles, as indicated in the figure. It is evi- 
dent from these results that the density ratio has the greatest 
effect upon the flow parameter due to the real gas effects. For 
a normal shock wave, y = 1.4 and infinite flow Mach number, 
the limiting density ratio is 6. In the present calculation with 
a free stream temperature of 810 R and a flow Mach number 
of 7.8, the density ratio across the normal shock wave is 10.25, 
as shown in Fig. 3. This larger density ratio causes the shock 
wave to move closer to the wedge surface as indicated in Fig. 
1, and causes the detached shock wave for a blunt body to 
move closer to the body surface (6). 


Experimental Procedure 


Hypersonic Shock Tunnel 


It is possible to check the calculations for the oblique shock 
wave for thermodynamic equilibrium, discussed in the previ- 
ous section, by using a two-dimensional wedge model in a hy- 
personic shock tunnel. A detailed description of this research 
equipment for producing hypersonic flows in the 24-in. diam 
test section at high stagnation temperatures is presented in 
(5). The free stream conditions in the test section for which 
the calculations were made correspond to a reservoir equi- 
librium stagnation temperature of 7200 R at a pressure of 500) 
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psi just ahead of the entrance to the hypersonic nozzle. Pre- 
viously it has been shown (6) that such a flow will be close to 
equilibrium conditions when expanded in the nozzle. It was 
for this reason that these particular upstream conditions were 
selected. 


Model and Instrumentation 


An adjustable, 60-deg included angle, wedge model with 
three pressure orifices on both surfaces was used in the in- 
vestigation. Wedge angles of 50 and 10 deg, 45 and 15 deg, 
40 and 20 deg, or 30 and 30 deg may be tested to obtain sur- 
face pressure for two different flow deflections at one time. 
The span of the model was 5 in. with the pressure orifices 
located near the center of the model. 

The Schlieren photograph (Fig. 4) of the shock wave was 
obtained with a single pass symmetrically arranged Schlieren 
system with a short duration condenser discharge spark for the 
light source. A delay circuit was used in conjunction with the 
Schlieren system to obtain Schlieren photographs at any de- 
sire time after the establishment of the flow. 

Piezoelectric barium-titanate gages were used to measure 
the surface pressure on the wedge. These were specifically 
constructed at the General Electric Research Laboratory for 
use in fundamental hypersonic research investigations. An 
orifice diameter of 7g in. was used to measure the pressure. 
This small diameter was adequate for the small barium- 
titanate gages. 

All the pressure gages were installed in the wedge model and 
then calibrated at the end of the shock tube portion of the 
hypersonic shock tunnel. A weak wave reflected from the 
end of the tube where the wedge model was installed was used 
to produce the increment in pressure for calibrating the gages. 
The shock Mach number for these calibration tests was about 
1.5, at which the agreement between the experimentally and 
analytically predicted pressures is quite good. It was found 
that the pressure gages must be dynamically calibrated in 
the model in order to obtain consistent and reliable pressure 
data in the hypersonic shock tunnel. 


Operating Procedure 


The model was then installed on the sting and positioned in 
the test section of the shock tunnel (see Fig. 4). The alumi- 
num diaphragm at the entrance to the conical nozzle was in- 
stalled, after which the dump tank housing the nozzle exit was 
evacuated to about 25 uw of mercury pressure. For a nozzle 
area ratio of 576, this dump tank pressure was low enough 
to minimize the strength of the starting shock wave. 

By adjusting the initial pressures in the driver and in the 
driven tube, the desired reflected pressures and stagnation 
temperatures at the end of the driven tube can be produced. 
The shock wave velocity is determined at 12 stations along 
the driven tube. At the end of the driven tube, a more ac- 
curate measurement of the shock velocity is made over a 
64-ft distance by means of a Berkeley counter. Thus, by de- 
termining the shock Mach number at the end of the tube and 
knowing the initial pressure and temperature of the air in the 
driven tube, the reflected equilibrium stagnation temperature 
can be calculated by using available thermodynamic informa- 
tion for air (7,8). It is assumed that equilibrium is attained 
rapidly after the shock wave. At high pressures the assump- 
tion of equilibrium after the reflected wave seems to be valid, 
based on recent results. The pressure after the reflected wave 
was measured with two pressure gages installed very close to 
the entrance of the conical nozzle. These reflected tempera- 
tures and pressures constitute the reservoir conditions for the 
convergent-divergent conical nozzle. 

Using these reservoir conditions, it is possible to calculate 
the pressure and temperature in the nozzle, assuming that the 
flow is always in equilibrium and isentropic during the expan- 
sion process. For a reflected pressure of 500 psi and an 
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Fig. 4 Schlieren photograph of adjustable wedge at 30 and 30 
deg angles from centerline 


equilibrium stagnation temperature of 7200 R the static and 
impact pressure along the centerline of the nozzle agreed very 
well with the value calculated for an equilibrium isentropic 
expansion (6). Thus, the initial free stream conditions as- 
sumed in the oblique shock wave calculations were actually the 
same as the results obtained for this nozzle. For the present 
investigation, the shock tunnel was operated with these re- 
servoir conditions to obtain the oblique shock wave angle and 
wedge surface pressure as functions of flow deflection. 


Results 


Shock Wave Angle 


For wedge angles less than 15 deg, measured from the free 
stream direction, the boundary layer becomes very thick be- 
cause of the low pressure downstream of the shock wave. In 
this case the outer edge of the boundary layer was displaced 
approximately 2.3 deg from the wedge surface. At lower 
wedge angles the flow deflection caused by the viscous effects 
becomes very large. For the wedge surface parallel to the 
free stream direction, the flow deflection may be as large as 
9 deg. This is due entirely to viscosity as determined from 
the shock wave experiments with a flat plate (9). A detailed 
discussion of these experiments along with the associated 
theory for the slip phenomena and the strong shock wave- 
boundary layer interaction is presented in (9 and 10). 

Since the flow deflection due to viscous effects is large for the 
range of 0- to 10-deg wedge deflection, no pressure information 
was obtained with the present wedge model. A flat plate 
model with pressure gages and with an adjustable angle of at- 
tack will be used to investigate the shock wave-boundary 
layer interaction at small flow deflections. In this case, the 
boundary layer displacement angle is the same order of mag- 
nitude as the surface deflection. To minimize these effects, 
the present wedge shock wave angles and surface pressures 
were obtained for wedge angles in the range 15 to 45 deg. The 
maximum angle was limited by the shock wave becoming 
curved as the flow detachment angle was approached for real 
gas conditions (Figs. 1 and 5). 
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Fig. 5 Comparison of experimental and theoretical shock wave 
angle as function of wedge angle 


The shock wave angles in Fig. 5 were determined from the 
Schlieren photographs (Fig. 4) and are plotted as a function 
of the wedge angle without any boundary layer correction. 
For all wedge angles the shock wave angles are greater than 
those calculated for the real gas case as shown in Fig. 5. If the 
wedge angle is corrected by the boundary layer displacement 
angle as estimated from the Schlieren results, the real gas 
theoretical shock wave angle approaches the observed shock 
wave angle as indicated in Fig. 5. The corrected results agree 
very well for the lower wedge angles, but become slightly 
lower at the higher wedge angles. It is evident by comparing 
the perfect gas case y = 1.4, presented in Fig. 1, with the ob- 
served shock angles, that the real gas effects decrease the shock 
angle for a given flow deflection. The angle for flow detach- 
ment is increased. Thus for the present test condition the 
assumption of flow being in thermodynamic equilibrium be- 
fore and downstream of the shock wave is a reasonable one 
based upon the shock wave angle corrected approximately for 
the boundary layer displacement. 

Further investigations will be conducted to obtain shock 
wave and pressure information in the flow deflection range 
from 45 to approximately 54 deg, using symmetrical wedge 
models. With the present wedge model, it was not possible to 
obtain reliable results in the vicinity of the shock wave de- 
tachment angle because of the unsymmetrical shape. For a 
flow Mach number of 7.8 at the test conditions, the flow de- 
flection angle for shock detachment was approximately 43.5 
deg for y = 1.4, and 54.5 deg for the real gas (Fig. 1). The 
observed shock wave angles in Fig. 5 indicate this delay in the 
detachment flow angle. 


Wedge Surface Pressure 


The experimental wedge surface pressure ratios were ob- 
tained in conjunction with the shock wave angles and under 
the same free stream condition. Corresponding to this free 
stream condition in the test section, the equilibrium stagna- 
tion temperature behind a normal shock wave was 6070 R. 
This stagnation temperature in the test section was lower than 
the 7200-R stagnation temperature at the entrance to the 
nozzle because of the shock wave in the real gas. The de- 
parture from a perfect gas in the test section was caused by 
the low pressure at high temperature. 

In Fig. 6 the wedge surface pressure ratio (surface pressure 
divided by the free stream pressure) for a wedge angle range 
of 15 to 45 deg is presented as a function of the wedge 
angle. The theoretical real gas pressure ratios were calculated 
by assuming the flow to be in equilibrium ahead and down- 
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stream of the oblique shock wave and using the thermo: y- 
namic data for air presented in (3 and 4). The caleulations 
were made assuming the flow to be inviscid; that is, the 
boundary layer on the wedge surface was neglected. The 
vertical bars for the pressure ratios indicate the approxim:te 
range of the pressure measured by the gages for a given wedge 
angle. 

The agreement between theoretical and experimental pres- 
sure results is quite satisfactory for the wedge angles con- 
sidered (Fig. 6). The agreement is approximately within the 
experimental scatter of the pressure gages. At a wedge angle 
of 43.5 deg, the shock wave will be detached for the case of 
perfect gas with y = 1.4, as indicated in Fig. 3. For this 
condition the pressures on the wedge will be close to the normal 
shock pressure ratio of 70. But, due to the real gas effects, 
the shock wave is attached at a 45-deg wedge deflection angle, 
and the surface pressure ratio is approximately 53. Thus for 
this moderately high stagnation temperature in the test sec- 
tion, the real gas effect upon the wedge surface pressure was 
appreciable for a large flow deflection. 


Conclusions 


It has been possible to measure the surface pressure for a 
wedge at several angles of attack for high stagnation tempera- 
ture by the use of barium-titanate pressure gages. 

The experimentally observed shock wave angles agreed 
with the theoretical values by assuming the flow to be in 
thermodynamic equilibrium ahead and downstream of the 
shock wave for a stagnation temperature of 6070 R after a 
normal shock wave and flow Mach number of 7.8. 

At low wedge angles the flow deflection caused by the 
boundary layer was appreciable, and, consequently, the ac- 
tual shock wave angle was much greater than the calculated 
inviscid value. 

The real gas effects delayed the flow deflection angle at 
which the shock wave becomes detached from 43.5 deg for a 
perfect gas with y = 1.4 to approximately 54.5 deg. 

The observed experimental surface pressure ratio agreed 
very closely with the theoretical results for equilibrium flow 
conditions. The pressure ratios at the large wedge deflection 
angles were appreciably less than that for the perfect gas 
y = 14 result. 
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Effect on a Rocket of the Oblateness 


of a Planet 


WILLIAM A. ALLEN! 


U. S. Naval Ordnance Test Station 
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An approximation technique is used to solve the mathematical problem of a rocket moving 
radially away from a planet. The acceleration of the rocket is given by the relation 0°r/0O(? + 
g(r) = c, where c is a constant and g(r) is the radial acceleration due to gravity at the distance r 
measured from the center of the planet. The planet is characterized by a gravitational potential 


that includes effects of its oblateness. 


A first integral, relating velocity and distance, is obtained 


immediately; a second integral, relating time and distance, is obtained by a quadrature. This 
second relation is expressible approximately in terms of elliptic integrals reducible to elementary 


types that have been tabulated. 


Inclusive of second-order terms the solution has been carried 


through completely, and the inclusion of fourth-order terms has been indicated. The approximate 
results are compared with the exact values obtained by numerical integration. Evaluation of a 
numerical example indicates that the solution inclusive of second-order terms is equivalent to the 
exact solution, and that the fourth- and other higher-order terms are negligible for the case of a 
rocket leaving Earth with an acceleration about 10 times that due to gravity. 


PAPER has appeared (1)? in which the problem was 
solved of a rocket that starts from rest, moves radially 
away from a planet and recedes to infinity against a force 
of gravity which diminishes as the inverse square of the dis- 
tance measured from the center of the planet. The motion 
of the rocket at a given distance r was determined by the 
requirement that its acceleration relative to the planet, plus 
the acceleration due to gravity at that distance, was equal 
to a constant value c. Thus, an accelerometer within the 
rocket would always record this constant value c; there would 
be no rate of change of radial acceleration within the rocket. 
The purpose of the present paper is to generalize the pre- 
vious results to the case where the external gravitational po- 
Received March 9, 1959. 
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tential of the planet can be characterized by the relation (2) 
3 5 
r 3 r 35 r 
[1] 


where P, = P,(cos 0) and Py = P;(cos 8) are the even Le- 
gendre polynomials defined by the relations 


P.(cos 8) = (1/2) (3 cos? 8 — 1) [2] 
P,(cos 8) = (1/8) (35 cost 8 — 30 cos? 8 + 3) [3] 


The quantity r in Equation [1] is distance measured from 
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the center of the planet; @ is the co-latitude; g is a constant 
that specifies the acceleration due to gravity at the distance 
a if the mass of the planet were concentrated at, or were 
arranged symmetrically about, its center. The radius of 
such a hypothetical planet must be less than or equal to a. 
In the case of an oblate planet, the quantity a is the equa- 
torial radius, and the two quantities J and D are coupling 
constants fitted to gravitational measurements on the surface 
of the planet. Equation [i] has the general form that ap- 
plies for a homogeneous oblate ellipsoid (3). Most planets 
can be considered, as a first approximation, to be homo- 
geneous oblate ellipsoids (4). Equation [1] has been used 
in recent papers to describe the effect of the oblateness of 
Earth on a near satellite (5,6,7). 


Mathematical Formulation and Solution 
of the Problem 


The problem reduces to the solution of the partial differ- 
ential equation 
r =c [4] 
The dots in Equation [4] signify differentiation with respect 
to time. Substitute Equation [1] into Equation [4], and 
make the substitutions s = r/a, k? = g/c, to obtain the 
acceleration of the rocket relative to the planet; this accelera- 
tion can be expressed by the relation 


§ = ca! {1 — k?[s~* — + (8/7) PwDs-*]} [5] 
Note that s =#/a and s=7/a. Integrate Equation [5], and 
introduce the boundary condition s = 0 when s = 1 to obtain 
the velocity of the rocket relative to the planet; this velocity 


Equation [7] is an improper integral. Integrate by parts to 
transform Equation [7] into the expression 


2a 1/2 s(s — 1) 

kt s(s — + 
2 J1 {s[s — k(1 + M)]}*” 


The quantity N in Equation [9] is specified by he relation 


N = — (2/83) PoJ (1 + 28-! + 387?) + 
(8/35) (1 + 28! + + 48-3 + [10] 


Equations [5, 6 and 9] constitute a rigorous formal solution 
of the problem. The integral of Equation [9], however, must 
be evaluated by numerical methods. 

A perturbation technique will now be applied to obtain an 
approximate analytical solution of Equation [7] The quan- 
tities J and D in Equation [8] are small for most planets in 
comparison with unity; therefore M is small in comparison 
with unity, and the quantity in square brackets in Equation 
[7] can be expanded by means of the binomial theorein. 
Equation [7] can be approximated by the relation 


[1 +(1/2)k2(s — k?)—1M +(1/2)(3/4)kXs — k?)-2M2 + .. |sds 
[4s(s — 1)(s — k*)]'” 


(11) 


Substitute the value of M from Equation [8] into Equation 
[11] and collect terms to obtain the relation 


8 (1 + s7! + 8-*)sds 


2a\'/2 s sds 1 4 
[4s(s — 1)(s — +m +4 ro) (s — k*)[48(s — 1)(s — + 


+ s—*)sds 


4 


35 1 (s — k*)[4s(s — 1)(s — k2)]'” 


can be expressed by the relation 
P 2c8 
— 1)(s — + 


3 “7 


Integrate Equation [6] to obtain the time elapsed for the 
rocket to reach the distance s; this interval can be expressed 
by the relation 


2a 


sin 


1 
+ 


s(1 + 2s-! + 38-2 + + s~*)sds 
(s — k*)?[4s(s — 1)(s — k?)] 


Equation [12] constitutes an approximate solution of Equa- 
tion [7] to terms of order D and J*. Equation [12] can be 
expanded into the sum of 11 integrals of the form 


I s™-" ds 

mn (8 — [48(s — 1)(8 — 
where m = 0,n = —1;m=1,n=0,1,...,4;m = 2,n = 1, 
2,...,5. By means of the substitution s = sin~? g, Equa- 
tion [13] can be transformed into the expression 


[13] 


sin?” gdy 
Inn = (1 — sin? [14] 


where g, = sin~! (s~'/?). Complete integrals of some of the 
forms indicated in Equation [14] have been tabulated (8). 
The incomplete integrals, however, must be reduced to simpler 
forms before evaluation can be effected. Integration of 
Equation [14] by parts yields the following recursion formula 


2n — 2 «/2 sin?"-? 


sin?" gdg 
(1 — k*sin? g)°/) (m+) (1 — k? sin? g)°/) 


2n —3 
k(2m —1) Je (1 — k*sin? k#(2m — 1) (s — k2)°/)@m—1) 


The quantity M in Equation [7] is specified by the relation 


M = — (2/3) (1 + + 8-*) + 
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~ K22m —1) (1 — sin? + 


an—4 1 m—n+'/2 1/2 
sin gde 8 (s ) [15] 


The integrals of Equation [14] can be evaluated by means of 
successive application of Equation [15] and the use of stand- 
ard reduction formulas (9). 

Consider the order of magnitude of the second- and fourth- 
order corrections to the zero-order solutions previously pub- 
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tute severe restrictions on the problem. In general, the 
corrections will be much smaller than those4ndicated. 


lished (1). The exact solutions, specified by Equations [5, 
6 and 9], can be written in the approximate forms 


(1 + =) [16] Evaluation to Second Order 
To 
In the case of the planet Earth, evaluation to fourth order 
inh ( ‘es Av % Atv ) (17] is not justified in view of experimental uncertainty in D and 
; 0 to J*. There is considerably less information about the effect of 
the oblateness of the other planets. Evaluation of Equation 
es (: + A*t x “*) [18] [12] will be carried out only to terms of J in order to illustrate 
: to ty the method. Further extension of the calculation to terms 
of order D and J? involve merely additional mathematical 
where %, % and t& are the zero-order solutions. The quan- manipulations of the same kind as those performed. Equa- 
a tion [12] can be written in terms of J in the form 
2a\'/2 dy | dy sin? gdg 
(3 ) sin? g(1 — k* sin? g)'/* 3 (1 — sin? ¢)’/? + (1 — k? sin? + 


sin‘ gdg 

The last integral of Equation [25] can be reduced by Equa- 
tion [15]; Equation [25] can then be written in the form 


2a\"/? dg 1 agg dg 
( ) \f sin? g(1 — k?sin? g)'/? at | (1 — sin? 


c 


(1 —k?sin? (1 — sin? + Je: (1 — k* sin? + k? \s(s — k?) [26] 


The elliptic integrals of Equation [17] have been tabulated 
(9,10). Equation [26] can be expressed in the form 


2a\"/? 1 2a\'/? ‘ 
t= (*) [K — F(¢1, k) — + k) + — — — 3 Pad (*") + kt) [K — F(g,k)] + — 


| '/2 
+ [BE — + 5 [= | + + 1)s — + 1) +.. [27] 


The respective quantities K and F(¢i,k) in Equation [27] 

tities A, A%, A% and A‘?, A‘, A are the respective are complete and incomplete elliptic integrals of the first 

second- and fourth-order differential corrections obtained by kind; the respective quantities Z and E(¢,,k) are the com- 

the methods used in this paper. The differential corrections plete and incomplete elliptic integrals of the second kind. 

will have their maximum effect relative to the zero-order The first term of Equation [26] is the result previously pub- 

effect when s = 1 and when c is comparable to g. Expand lished (1); the remaining terms constitute the new results 
derived in this paper. 


Equations [5, 6 and 9], and evaluate at the point s = 1 to ob- 
tain the relations 
(c [: 8 ( 9 ) (19) Numerical Example 
A numerical example has been evaluated in order to com- 
: 2a(s — 1)(cs — g) |'? g pare the approximate and the exact solutions. Recall that 
a Pt the boundary conditions require that the rocket start from 
if % 1 ¢ \ rest at the equatorial distance a from the center of the planet. 
= ( —£) P,D += () PJ ‘| [20] This requirement means that the rocket can start from rest 
. = from the surface of the planet at the Equator only. At 
—1)]” 9 \py4 other latitudes (see Fig. 1) the rocket must start from an 
cs —g elevation just sufficient to compensate for the difference in 
ass 3f¢\3 distance to the center of the planet caused by its oblateness. 
3 ( ) PD + - (-) PJ ‘| [21] Launching from considerable altitude by means of a plane or 
—— 2\e~9 a balloon is definitely an engineering possibility. There are 
Set Po = 1, Ps = 1 and c = 2g in Equations [19, 20 and 21], considerable advantages to an air launched rocket, since 
and use the values of J and D from Table 1 to obtain the about 94 per cent of the atmosphere lies below 20,000-m 
following relations that apply to Earth altitude (11). 
7 = [e — (g/s*)][1 + 0.003274 — 0.000012] [22] Launching at a latitude other than 2/2 requires that the 
; [ 8 ] ciniintnaiiities » trajectory in order to cancel the effect of rotation. At lati- 
2Qas(s — 1) ]'/2 tudes other than zero and 7/2, the rocket must also undergo 
t= esa [1 — 0.001637 + 0.000010] [24] a slight acceleration normal to its trajectory in order to re- 
OB: main on a radius vector, since the force of gravitation given 


Equations [22, 23 and 24] can be used to estimate the magni- by VU is not central. Reference to Equation [1] suggests 
tude of the second- and fourth-order effects of oblateness. that the maximum perturbation effect occurs for a rocket 
The values used to obtain Equations [22, 23 and 24] consti- launched at a pole where the effect is twice that, and in oppo- 
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site direction to, the effect produced at the Equator. Refer- 
ence to Equation [27] suggests that the second-order effects of 
the oblateness of a planet vanish entirely at the latitude 
35°-15'-52”. This result follows from the fact that P2(cos @) 
has a root at this value. Similarly, the fourth-order ef- 
fects vanish at latitudes 19°-52’-33” and 59°-26’-40". 

In view of the preceding discussion, the problem has been 
evaluated for a rocket launched at the altitude a — b = 
21,477 m above a pole of Earth; see Fig. 1. Values of the 


ROCKET 


START 
PLANET— 


Fig. 1 Rocket moving radially away from an oblate planet. 
Dashed circle indicates boundary conditions 


parameters used in the solution are listed in Table 1. The 
constant g in Table 1 was determined so that the acceleration 
due to gravity at the pole is given by Equation [1]; that is, 
g{(a/b)? — 2J(a/b)* + (8/7)D(a/b)*] is the theoretical 
value of gravity at the pole (983.2213 cm per sec?). The 
values of the equatorial and polar radii of Earth were taken 
to be a = 6,378,388 m and b = 6,356,911 m, respectively. 
The polar axis b is a computed quantity based upon the as- 
sumption that a/(a — 6) = 297.0, exactly. If Earth were a 
homogeneous ellipsoid, the coupling constants J and D would 
be functions of the eccentricity e of Earth; that is, J would be 
(3/10)e? and D would be (3/8)e* where e? = 1 — (6/a)’. 
The theoretical value of the coupling constants for a homo- 
geneous ellipsoid are: J = 2.017 X 10-* and D = 16.9 x 
10-*. The experimental values J = (1.637 + 0.004) x 10°* 
and D = 10.6 X 10~* suggest that the actual gravitation: | 
field of Earth is somewhere intermediate between that of a 
homogeneous ellipsoid and a sphere. Since the mean density 
of Earth (5.517 g per cm‘) is considerably greater than the 
average density of the continents (2.679 g per cm*), Eart!: 
cannot be homogeneous. The value chosen for c in the last 
column of Table 1 was about 10.5 g. This magnitude of 
acceleration, chosen for mathematical convenience, represents 
a reasonable upper limit to the intensity of acceleration that 
can be sustained for a short time by a well-conditioned humai: 
being. 

Table 2 contains the results of the zero-order solution ob- 
tained from Equations [5, 6 and 9] by setting J = 0 and D = 
0, and using the value of the other parameters listed in Table 
1. The exact values, also listed in Table 2, were obtained 
from Equations [5, 6 and 9] by the use of all of the parameters 
listed in Table 1. The difference between an entry for the 
zero-order and the corresponding entry for the exact value is 
the sum of differential corrections of the second- and higher 
even orders. These total differential corrections are listed 
in Table 3 and are plotted in Fig. 2. The total differential 
corrections do not differ, to the precision of the numbers 
tabulated, from the second-order differential corrections. 
Thus, in this example, the second-order solution is equivalent 
to the exact solution. 


3 See (11), pp. 2-100. The geodetic quantities chosen in this 
paper are those of the International Ellipsoid adopted in 1930 by 
the International Geodetic and Geophysical Union. 


Gravity Radius Co-latitude 
g, m/sec? a,m 6, deg 
9.798, 288 6,378,388 0 


Table 1 Parameters used in the problem of a rocket launched from a pole of Earth 


Acceleration 
c, m/sec? 
102.609 


Coupling constants 


1.637 X 1073 10.6 107% 


Distance s, Velocity 7, km/sec 


Earth radii zero-order exact 
1.0 00.000 00.000 
10.933 10.934 
1.2 15.522 15.524 
1.3 19.074 19.076 
1.4 22.087 22.089 
1.5 24.755 24.756 
1.6 27.175 27.177 
29.409 29.409 
1.8 31.489 31.491 
1.9 33.449 33.450 
2.0 35.305 35.306 


Table 2 Zero-order and exact values for velocity, acceleration and elapsed time of a typical rocket 


Acceleration 7, m/sec? Time ¢, sec 

zero-order exact zero-order exact 
92.810 92.842 000.000 000.000 
94.511 94.533 117.044 117.025 
95.804 95.820 165.282 165.257 
96.811 96. 822 202.160 202.131 
97.609 97.618 233.154 233.122 
98 . 254 98. 260 260.389 260.355 
98.781 98.786 284.955 284.919 
99.218 99 . 222 307 .501 307 . 463 
99 . 584 99 .587 328.449 328.410 
99 . 894 99 .897 348 .094 348.054 

100. 159 100.161 366.648 366 . 607 
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File 3 Differential corrections for velocity, acceleration 
and elapsed time of a typical rocket - 
Distance s, Velocity Acceleration Time 
Earth correction correction correction 
radii Av, m/sec Ar, m/sec? At, sec 
1.0 0.00 0.031 —0.000 
1.2 1.84 0.015 —0.025 
1.4 1.95 0.008 —0.032 
1.6 1.89 0.004 —0.036 
1.8 1.79 0.003 —0.038 
2.0 1.68 0.002 —0.040 
2.2 1.59 0.001 —0.042 
2.4 1.50 0.001 —0.043 
2.6 1.42 0.000 —0.044 
2.8 1.36 0.000 —0.045 
3.0 1.30 0.000 —0.046 
Conclusions 


he differential corrrection that pertains to the J coupling 
constant is relatively small for the illustrative example 
studied, and the differential correction that pertains to the 
D coupling constant is negligible. The respective accelera- 
tion and velocity differential corrections Ar and A? approach 
zero as 8 approaches infinity. Fig. 2 indicates that the ve- 
locity correction A? attains a maximum value r = 1.958 m 
per sec in the neighborhood of the distance s = 1.37. The 
time correction At approaches a maximum value t = —0.05953 
sec as s approaches infinity. 

The restricted problem solved in this paper is characterized 
by an analytical solution. This result presumably can be 
generalized by perturbation techniques to the more realistic 
case where the acceleration and trajectory may vary from 
the assumed conditions. Results of the restricted problem, 
however, suggest that the corrections due to the oblateness 
of a planet are generally small or negligible. This conclusion 
can be applied directly to the more general case. 
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Mass Transfer and 
Shock Generated Vorticity 


H. HOSHIZAKT' 


Lockheed Aircraft Corp. 
Sunnyvale, Calif. 


The effects of shock generated vorticity and mass transfer on stagnation point heat transfer raies 
are investigated. The complete incompressible Navier-Stokes equations are considered in the flow 
region between the bow shock and the surface of spheres and cylinders. Boundary conditions are 
applied immediately behind the shock and at the wall. The numerical solutions to the flow equa- 
tions with air injected into the shock layer, show that the interaction between the vorticity gencr- 
ated by the wall and by the curved shock reduces the effectiveness of mass transfer cooling. The 
reduction in heat transfer rates due to mass injection is substantially less than predicted by bound- 
ary layer theory at Reynolds numbers below 10° for spheres and 10‘ for cylinders. The heat transier 
rates with mass injection were found to be 200 to 300 per cent greater than corresponding boundary 
layer values for the extreme cases investigated (strong shock, Re = 10° and large air injection rate;). 
These heat transfer rates are, however, less than the zero mass injection values. It is also shown for 
a sphere and for Pr = 1.0 that the increase in heat transfer is primarily dependent on the inverse of 


the difference between the vorticity at the wall and the vorticity immediately behind the shock. 


ASS transfer is important in several methods for pro- 
tecting the surface of hypervelocity vehicles. For ex- 
ample, on an ablating body the interface temperature of the 
liquid sublayer may be sufficiently high to vaporize the liquid. 
The vaporized mass will subsequently be injected into the air 
boundary layer. It has been shown by numerous authors 
(1-5)? that the injection of mass into boundary layers results 
in substantial reductions in conductive heat transfer rates. 
Thus, it is necessary to know the quantitative reduction in 
heat transfer rates with mass injection in order to accurately 
compute ablation rates. For a more detailed discussion of 
ablation the reader is referred to (6 and 7). 

Another method of protecting hypervelocity vehicles from 
severe aerodynamic heating is to use materials, such as graph- 
ite and Teflon which sublime at high temperatures. The 
sublimated mass is injected into the boundary layer, and again 
one must know the effects of mass injection on heat transfer 
rates to compute sublimation rates. 

Mass transfer cooling, the process whereby fluid is injected 
into the boundary layer from a porous wall, is also an effective 
means of alleviating aerodynamic heating. Obviously for 
such a scheme one must know the effects of mass injection on 
heat transfer rates. 

It is the purpose of this paper to show that the reduction 
in stagnation point heat rates by mass injection is markedly 
affected by shock generated vorticity for Re ~ 10° and lower. 
Reference (8) contains an analysis of the effect of shock 
generated vorticity on stagnation point heat transfer rates on 
a hemisphere, without mass injection, at low Reynolds num- 
bers. The results show that the interaction between the vor- 
ticity generated by the bow shock and the vorticity at the 
wall substantially increases the heat transfer rates over the 
values predicted by classical boundary layer theory. It was 
also shown for the special case of a sphere and for Pr = 1.0, that 
the increase in heat transfer rates is primarily dependent on 
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the inverse of the difference between the vorticity at the wall 
¢. and the vorticity immediately behind the shock ¢,. At 
high Reynolds numbers (Re > 10°), ¢. is much greater than 
¢, and the heat transfer rates are little affected by shock 
generated vorticity. As the Reynolds number is decreased, 
¢» decreases whereas ¢, is only slightly affected. Hence, at 
low Reynolds numbers, the heat transfer rates are increased 
over boundary layer values by virtue of the decreased mag- 
nitude of the quantity —¢,). Since boundary layer 
solutions show that mass injection decreases ¢, [see (5) for 
example], one would expect heat transfer rates considerably 
higher than those predicted by boundary layer theory when 
mass is injected into the shock layer at moderately low Rey- 
nolds numbers. 

These qualitative conclusions are verified by the results of 
the present analysis. Numerical solutions to the incompres- 
sible Navier-Stokes equations with air injection into the shock 
layer were obtained for flows in the vicinity of the stagnation 
point of hemispheres and cylinders. These solutions give 
heat transfer rates which are as much as 200 to 300 per cent 
higher than those predicted by boundary layer theory for the 
same injection rates. Although the present solutions are for 
single component flows, one can logically expect similar in- 
creases in heat transfer rates over the corresponding boundary 
layer solutions for multicomponent flows, which are more ap- 
propriate for ablation and sublimation calculations, 


Analysis 
Basic Equations and Boundary Conditions 


The present analysis is confined to the vicinity of the for- 
ward stagnation points of hemispheres and cylinders where 
the flow is assumed to be incompressible. The results of an 
incompressible analysis when referenced to the proper bound- 
ary layer results will give the correct trends when applied to 
a compressible flow. The magnitudes of the results of such 
an analysis should also be very good, since the difference in 
compressibility effects between a shock layer and boundary 
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layer analysis at moderately low Reynolds number is ex- 
pected to be a second-order effect. 

In this analysis the region between the body and bow shock 
(shock layer) is considered to be completely viscous, whereas 
the shock wave is assumed to be a discontinuity. Adams and 
Probstein (9) have shown by an order-of-magnitude analysis 
that these conditions can be expected if the ratio of body 
radius to the mean free path at the wall r./Aw, which is pro- 
portional to M.,Re,, is greater than 0(10)*. It was also 
shown for M,,Re,, > 0(10)? that the usual vorticity trans- 
port, continuity and energy equations are applicable. In this 
paper, the development of these equations closely follows the 
analysis of (8). However, for completeness, pertinent por- 
tions of the analysis will be repeated here with a slight gen- 
eralization to include stagnation points on cylinders aligned 
normal to the flow. It should also be noted that the analysis 
can be extended to include other bodies of the conic family. 

The vorticity transport, continuity and energy equations 
in the coordinate system shown in Fig. 1 can be written as: 

Vorticity transport 


(£) + v sin = 


Continuity 


fe) 
n+ in” = ¢ 
+1 sin" @u) + 30 (r™ sin” v) = 0 [2] 


Energy 


or . vor 1 oO oT 
(« or r E or (- 7) 
1 oT 
r? sin" 00 (sin | 


where n = 1 for the sphere, n = 0 for the cylinder, and the 
vorticity ¢ is defined by 


1] 0 ou 
r= 2/2 [4] 


The dissipation and pressure gradient terms in the energy 
equation have been omitted, since these terms were shown in 
(8) to be negligible. 

A stream function is now introduced and is assumed to have 
the form 


= (Usre" +3/V Re) +! Oh(n) [5] 
where 
= [(r/re) 1]V Re [6] 


In order to satisfy continuity, the velocity components are 
written as 


v = (U, sin 0/7") h(n) [7] 
u = —[(n + 1)U., cos 6/9" + Re] h(n) [8] 

where 
1 =1/rTw [9] 


and the prime denotes differentiation with respect to 7. 
Furthermore, let 
Ty” +1 é 


Kn) = = 
V Re sin 


1 h” 


(l—n)h'’ (n+ 1)h 
Re Re | [10] 


1960 


Fig. 1 Coordinate system 


gn) = (T Tw)/(Ts — Tw) [11] 


By means of these relations the vorticity and energy equa- 
tions become: 
Vorticity 


—(n + + (1 — n)h'= = 
yn +.) = 
(i + V Re, (1 + n)n Re 


Energy 
—(n + = [13] 


Although the vorticity equation is a fourth-order equation 
five boundary conditions are necessary to specify a solution 
completely. The additional boundary condition is used to 
determine the unknown shock detachment distance. Two ad- 
ditional boundary conditions on the temperature profile are 
necessary to solve the energy equation. These boundary 
conditions are given at the body surface or immediately be- 
hind the bow shock, which is assumed to be a concentric 
sphere or cylinder. At the wall 


[14a] 
[14b] 
[14c] 


and immediately behind the shock 


v, = U, sin 8 {14d] 
us = —pU,, cos [14e] 
T=T, [14f] 


where the strong shock approximations were used to obtain 
Equations [14d and l4e]. The @ momentum equation 
evaluated immediately behind the shock provides the addi- 
tional boundary condition necessary to determine the shock 
detachment distance. The @ momentum equation is 


vo O(r v) ou 
where 
ps = (1 — p)ppU.,” cos? [16] 
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Transformation of these boundary conditions results in 


h'(0) = 0 [17a] 

g(0) = 0 [17¢] 

h'(ms) = 7." [17d] 

h(n.) Re/(n + 1) [17e] 
gms) = 1 [17f] 


The @ momentum equation then becomes 


(n-—1)- 


Equation [25] it is seen that this reduction in ¢, and 0¢,,/dr 
tends to reduce the heat transfer rate, but there is a compen- 
sating effect since ¢, also appears in the denominator. The 
compensating effect increases as ruf. —> r.£.. This can be 
seen more clearly by rewriting Equation [25] in terms of /)(7) 
and the Nusselt number, and taking the limit as Re—> 


Nu 2[—h'"(0) + 2h’(0)/W Re — 8h(0)/Re*!?] 


VRe (h'(0) — 2h(0)/Re] — — 2h(m.)/Re] 


27.2" 


= [ Re | [2501 + Re (1 — p) — (n — | [18] 


Heat Transfer Relations 


The temperature distribution through the shock layer is 
obtained by integrating the energy equation. Integration 
of Equation [13] results in 


Prh 1 
g(n) = % exp | - f, (n + 1) + 
[19] 


where the constant of integration is 


” Prh 1 


[20] 
The heat transfer rate at the wall is given by 
[21a] 
Gu = Re/re)aT, — Tw) [21b) 
or in nondimensional form 
Nu/V Re = 2a [22] 
where 
Nu = hD/k [23] 
h = qu/(T; — Tw) [24] 


In the most general case, the constant of integration ap can 
only be evaluated after a solution of the vorticity equation 
has been obtained. However, for the special case of a sphere 
and Pr = 1, Equation [20] can be integrated to obtain the 
heat transfer rate in terms of quantities evaluated at the wall 
and immediately behind the shock wave. The integration 
for the zero mass injection case is presented in (8). In order 
to include the effects of air injection, the wall boundary con- 
dition on the u velocity component is changed from u = 0 to 
u = Uy. For the heat transfer rate at the wall there results 


K( —Tw o¢../or + — Tw) 
rss) + sin + Pp) + Uy tan 6 


Gu [25] 
The effect of shock generated vorticity on the heat transfer 
rate is clearly shown by Equation [25]. At large Reynolds 
numbers, rf» is much greater than r.{,, so that the vorticity 
behind the shock has little effect on the heat transfer rate. 
At low Reynolds numbers, O(r..) is equal to O(r.¢.), and 
the heat transfer rate will differ markedly depending on 
whether or not shock generated vorticity is taken into ac- 
count. 

The main effect of mass transfer is in reducing the vorticity 
({~) and the vorticity gradient (0f,/0r) at the wall. From 
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Equation [26] shows that the direct effect of the mass inj:c- 
tion velocity [u» a h(0)] is relatively small. The limit of 
Equation [26] as Re— o is 


Nu —2h'""(0) 
V Re h''(0) h''(ns) 


which can be shown to be equivalent to the result obtained 
from boundary layer theory when applied to flow fields with 
external vorticity. Now, h’’’(0) is proportional to 0f¢,,/0r, 
h''(O) roftw and « refs, and since mass injection 
creases h’’’(0) and h’’(0) but has little effect on h’’(n,), it is 
easily seen that the effectiveness of mass injection in reducing 
heat transfer decreases as h’’(0) — h’’(n,) (or as the Reynolds 
number decreases). One other observation can be made about 
the effect of mass injection on heat transfer. Since the increase 
in heat transfer rates resulting from shock generated vorticity 
can essentially be attributed to the term (ruf. — rf.) in 
Equation [25], and since mass injection decreases ¢,, one 
would expect that the effect of shock generated vorticity on 
heat transfer rates at a given Reynolds number is larger with 
injection than without injection. 


Numerical Results 


Numerical solutions to the vorticity transport equation 
(Eq. [12]) were obtained by successive approximations. 
Equation [12] was first put into integral form by repeated 
integration. The constants of integration were then evalu- 
ated by means of the boundary conditions on the velocities 
(Kq. [17]). 

Estimates were first made of the shock detachment dis- 
tance 7, and the function h(n) for specified values of the 
density ratio p, the Reynolds number Re, and the injection 
rate parameter h(0). The integral form of the vorticity equa- 
tion was then solved, holding 7, constant, by successive ap- 
proximations. The function h(n) and its first three deriva- 
tives were required to converge to within 0.001. The solution 
thus obtained was then used to evaluate the error in the 6 
momentum equation boundary condition (Eq. [18]). This 
procedure was repeated for various values of 7, until Equation 
[18] was satisfied. 

Representative numerical results for the sphere and cylinder 
are presented in this section. Shock layer velocity, tempera- 
ture and vorticity profiles along with various Nusselt number 
ratios are discussed. Additional numerical results for the 
sphere and cylinder are presented in (10). 


Shock Layer Profiles and Shock Detachment Distance 


Velocity, temperature and vorticity profiles in the shock 
layer are shown in Fig. 2 for both the sphere and the cylinder. 
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Fig. 2 Shock layer profiles for | = 0.06 


The sphere and cylinder velocity profiles away from the wall 
resemble the inviscid profile obtained by Lighthill (11) for 
the sphere and by Whitham (12) for the cylinder. The injec- 
tion of mass is seen to reduce the velocity gradient at the wall 
and to thicken the shock layer. The temperature profiles 
show that mass injection reduces the temperature gradient 
at the wall, and the vorticity profiles show that mass injection 
reduces the vorticity at the wall. The vorticity immediately 
behind the bow shock is affected little by mass injection or by 
a reduction in Reynolds number. 

Shock detachment distances are shown in Fig. 3 as a func- 
tion of the Reynolds number and the injection rate parameter, 
for p = 0.06. The shock detachment distance is seen to in- 
crease with decreasing Reynolds number. The injection of 
mass into the shock layer also increases the shock detachment 
distance. 


Nusselt Number Ratios 


All of the heat transfer data are presented in terms of Nus- 
selt number ratios. The boundary layer Nusselt number, 
which is used as a reference value, was obtained from (11 and 
13) for the zero mass injection case for the sphere. For the 
sphere mass injection cases, the appropriate boundary layer 
equations were solved using the digital computer program of 
(5). The boundary layer Nusselt numbers for the cylinder 
were obtained from (12 and 14). 

Ratios of shock layer theory Nusselt numbers, to the zero 
mass injection boundary layer theory Nusselt numbers, are 
presented in Fig. 4 for both the sphere and cylinder, for p = 
0.06. The data are presented for various Reynolds numbers 
and are plotted as a function of the injection rate parameter 
h(0). These curves show the reduction in Nusselt number 
due to mass injection as predicted by shock layer theory. 
The large effect of vorticity on the sphere with zero mass in- 
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Fig. 3 Shock detachment distances 


jection causes the Nusselt numbers at low Reynolds numbers 
and small injection rates to be higher than the zero mass in- 
jection boundary layer Nusselt number. The increase in 
Nusselt numbers with zero mass injection accounts for a major 
portion of the percentage difference between shock layer 
theory and boundary layer theory. The increase in Nusselt 
numbers owing to shock generated vorticity is much smaller 
for the cylinder than for the sphere. 

The reduction in Nusselt numbers owing to mass injection 
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Fig. 4 Ratio of shock layer theory Nusselt number to zero injection boundary layer theory Nusselt number 


from the zero mass injection shock layer values is shown in 
Fig. 5. These curves are similar to those obtained from 
boundary layer theory. The difference between the curves 
obtained from shock layer theory and those obtained from 
boundary layer theory gives the increased effect of shock gen- 
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Fig.5 Nusselt number reduction as a function of mass injection 
tate according to shock layer theory, 6 = 0.06 
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erated vorticity when mass transfer is present. This can be 
seen more easily by considering the ratio of shock layer to 
boundary layer Nusselt numbers for an arbitrary value of the 
injection rate parameter. This ratio can be written as 


(Nu)sz (Nu)s1,0 [ | 
(Nu)sr L(Nu)ar/(Nu) 


The bracketed quantity is equal to the ratio of the shock layer 
to the boundary layer curve given in Fig. 5. As an example, 
for Re = 10%, p = 0.06 and h(0) = —0.50, this ratio for the 
sphere is equal to 1.30. In other words, under these condi- 
tions, the effect of shock generated vorticity on heat transfer 
is increased 30 per cent by mass injection. Again for the 
sphere and the same conditions, the ratio of shock layer Nus- 
selt number to boundary layer Nusselt number for zero mass 
injection is 1.63 (Fig. 4). Therefore, for an injection rate 
parameter of —0.50, shock layer theory gives a Nusselt num- 
ber at the stagnation point of a sphere which is 2.12 times as 
large as the Nusselt number predicted by boundary layer 
theory for the corresponding value of the injection rate pa- 
rameter. Note that for a cylinder, the increased effect of 
shock generated vorticity when mass transfer is present is 
less than it is for a sphere. 

The increase in Nusselt number over the corresponding 
boundary layer value as predicted by shock layer theory is 
shown in Fig. 6 as a function of Reynolds number for various 
injection rates. The increase in Nusselt number becomes 
significant for Reynolds numbers less than 105 for spheres and 
for Reynolds numbers less than 10‘ for cylinders. At low 
Reynolds numbers and high mass injection rates, shock layer 
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theory for a sphere predicts Nusselt numbers which are about 
twice the values predicted by boundary layer theory for the 
corresponding conditions. 


Application of Results 


The results of the present analysis were used to calculate 
the heating history of a vehicle re-entering from a satellite 
orbit. The trajectory used in these calculations is shown in 
Fig. 7. This trajectory is for W/CyA = 100 lb/ft? and L/D 
= 0.50. The initial re-entry angle is — 1.0 deg, and the initial 
velocity is 25,400 fps. Note that the vehicle begins to skip at 
about 400 sec. 

The heat rates at the stagnation point of a 1.0-ft diameter 
hemisphere, based on the trajectory shown in Fig. 7, are pre- 
sented in Fig. 8 for various injection rates. Very early in the 
trajectory, shock layer theory predicts heat rates which are 
higher than the free molecule values. This indicates that the 
present shock layer theory is not valid at these low Reynolds 
numbers. In Fig. 8, the shock layer and boundary layer heat 
rates have been extrapolated to the free molecule values to 
obtain a smooth transition from continuum to free molecule 
flow. 

Karly in the trajectory, when the Reynolds numbers are 
low, the heat rates predicted by shock layer theory are ap- 
proximately 50 per cent higher than the values predicted by 
boundary layer theory. Very late in the trajectory, the heat 
rates predicted by shock layer and boundary layer theory be- 
come equal. The total heat input predicted by shock layer 
theory for an injection rate parameter equal to —0.50 is 24 
per cent greater than the total heat input predicted by bound- 
ary layer theory. It should be noted that an injection rate 
parameter equal to —0.50 corresponds to an injection rate of 
a few hundred Ilb/sec-ft? under satellite re-entry conditions. 


Nomenclature 

C, = Specific heat at constant pressure 

D = body diameter 

f = modified stream function 

g = nondimensional temperature, (77 — 7'»)/(7's — Tw) 

h = modified stream function, also heat transfer coefficient 
k = thermal conductivity 

M = Mach number 
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coordinate normal to body 


Tr = 
Re = Reynolds number, U.rw/vs 

T» = body radius 

rs = shock radius 

s = nondimensional temperature 

T = temperature 

u = velocity component in r direction 

U.. = free stream velocity 

v = velocity component in 6 direction 

6 = shock detachment distance 

¢ = vorticity component 

7 = nondimensional normal distance 

= 

@ = spherical coordinate, body angle 

> = mean free path 

uw = dynamic viscosity 

vy = kinematic viscosity 

— = (see Equation [10]) 

p = density 

p = density ratio across normal shock, p./ps 
¢ = spherical coordinate (circumferential angle) 
y = stream function 

Subscripts 

BL = boundary layer 

QO = zero mass injection 

s = conditions behind shock wave 

SL = shock layer 

w = conditions at body surface 

o = free stream conditions 


Primes denote differentiation. 
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Technical Notes 


Thrust Orientation Patterns for Orbit 
Adjustment of Low Thrust Vehicles 


H. BROWN! and J. R. NELSON? 


Flight Propulsion Laboratory, General Electric Co., 
Evendale, Ohio 


Low thrust propulsion systems can be used for the pre- 
cision adjustment of space vehicle orbit characteristics. 
The most effective application of these low levels of thrust 
requires the use of thrust orientation patterns peculiar to 
each type of mission. The thrust patterns most suitable 
for altitude correction, eccentricity reduction and Earth 
escupe missions are identified and typical propulsion times 
indicated for each. 


GPAck vehicles powered by low thrust propulsion systems 
—heat transfer, plasma and ion rockets—operate in a 
gravitational field orders of magnitude greater than the applied 
thrust. The path of the vehicle will at first closely resemble 
the conic section of an unperturbed orbit. The elements of 
the conic section, however, will be continuously changing due 
to the action of the applied thrust. Propulsion effects can be 
segregated from gravitational effects most effectively by a 
direct numerical integration of the osculating orbit elements 
(1).3 

Orbit Characteristics 


Six elements are required to completely define an unper- 
turbed orbit in inertial space. A typical set of elements con- 
sists of: The semi-latus rectum and the eccentricity which 
define the geometry of the orbit, the inclination angle and the 
longitude of the node which define the orientation of the orbit 
plane, and the argument of the perigee and the time of perigee 
passage which define the orientation of the orbit in its plane. 
The location of the vehicle in its orbit is given by the true 
anomaly @ which is a function of the orbit elements and the 
time. The thrust orientation angle \ will be measured from 
the radius vector to the thrust vector in the direction of rota- 
tion. It will be assumed to lie in the orbit plane. 

The orbit elements are constant unless propulsive thrust or 
other perturbations are present. Except for short period 
oscillations, the elements defining the geometry of the orbit 
will be effected only by the propulsive thrust. 


Resonance 
The ability of the thrust to produce desired changes in the 
orbit elements is dependent predominantly upon the true 
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anomaly and the thrust orientation angle and also upon the 
values of the elements and the thrust to local vehicle weight 
ratio. The rates of change of each of the elements can be 
represented in the A-¢ plane as illustrated in Fig. 1 for the 
semi-latus rectum and the eccentricity. Note that the maxi- 
mum rates of change occur at the midpoints of the positive 
and negative regions. 

It is evident that the use of transverse thrust \ = 90 deg 
will result in a continuous increase in the semi-latus rectum 
and in oscillations in the eccentricity. Similar oscillations 
will occur in the argument of the perigee. Both the amplitude 
and the period of these oscillations will increase as the semi- 
latus rectum increases. A stable “resonance” condition is 
ultimately reached when the rate of change of the argument of 
the perigee is equal to the angular velocity of the vehicle. 
At this point, the true anomaly remains essentially constant— 
at a value less than 90 deg — and the eccentricity increases 
continuously. Similar resonance phenomena have been en- 
countered with many of the other thrust patterns investigated. 


Thrust Orientation Patterns 


The orbit element characteristics illustrated in Fig. 1 
can be used to define specific \-¢ relationships which will 
maximize the rate of change of a particular orbit characteris- 
tic. The A-@ relationship for some of the thrust orientation 
patterns studied are illustrated in Fig. 2. These are: 

1 Tangential, maximizes the rate of increase of total 
energy (2). 

2 Transverse, maximizes the rate of increase of semi-latus 
rectum, a measure of angular momentum (3). 

3 Radial, dissipates excess kinetic energy thereby provid- 
ing a reduction in eccentricity (4,5). 

4 Circularizing, maximizes the rate of decrease of orbit 
eccentricity. 

5 Constant Eccentricity, provides a continual increase in 
semi-latus rectum without the increase in eccentricity that 
accompanies the use of the transverse thrust pattern. 

6 Constant Perigee Altitude, provides a continual increase 
in orbit energy without the increase in perigee potential energy 
that accompanies the use of the tangential thrust pattern. 

7 Maximum Eccentricity, the opposite of the circularizing 
thrust pattern. 

The precise shape of the curves is dependent upon the 
physical shape of the orbit as indicated by the instantaneous 
eccentricity. 


Earth Escape Mission Performance 


The Earth escape mission shall be assumed to consist of 
the alteration of an initial low altitude-circular orbit to a para- 
bolic orbit. The pattern identified for this mission will be 
suitable for both lunar and interplanetary missions. The 
propulsion system is required to increase both the eccentricity 
and the total energy of the orbit to the parabolic level. 
Atmospheric drag must be minimized by maintaining a con- 
stant or increasing perigee altitude. 

Fig. 3 contains the results of escape mission calculations 


Eprror’s Nore: 


The Technical Notes and Technical Comments sections of ARS JourNAL are open to short manuscripts describing 


new developments or offering comments on papers previously published. Such manuscripts are usually published without editorial 


review within a few months of the date of receipt. 


JuLty 1960 


Requirements as to style are the same as for regular contributions (see masthead page). 


635 


SEMI-LATUS RECTUM,P 


° 


T1180 


© 


THRUST ORIENTATION ANGLE, 


Oo 90 180 270 360 
TRUE ANOMALY, 


ECCENTRICITY, e 


ti80 


THRUST ORIENTATION ANGLE, 


fe) 90 180 270 360 
TRUE ANOMALY, 


g.1 Orbit element variation, in thrust orientation angle true 
anomaly plane 


Fi 


CIRCULARIZING 


t180 


CONSTANT’> 
ECCENTRICITY 


20 


CONSTANT 
IGEE 
PERIGE 


Fig. 2 Thrust orientation programs, in thrust orientation angle 
true anomaly plane 


636 


for a thrust to weight ratio of (10)~*. The characteristic 
variation of eccentricity with semi-latus rectum is illustrated 
for the transverse, tangential, maximum eccentricity and 
constant perigee thrust orientation patterns. The most 
effective pattern for achieving the escape condition appears to 
be tangential thrust orientation. 

A propulsion period of 12.25 hr is required for achieving an 
eccentricity of 1 from an initial 300-mile orbit with tangen- 
tial thrust and a thrust to weight ratio of (10)~%. Calcula- 
tions for thrust to weight ratios of (10)—* and (10) ~ indic:ite 
corresponding times of 153 and 1886 hr, respectively. Escape 
appears to occur at the semi-latus rectum at which the thrust 
to local vehicle weight ratio is approximately 1 at 35,000, 
110,000 and 375,000 miles for (10)~?, (10)-* and (10)~4, 
respectively. This is due to the resonance requirements. 

Although the eccentricity varies continuously during the 
mission, the actual path of the vehicle closely approximates a 
circular spiral. 


Eccentricity Correction Mission 


The eccentricity correction mission shall be assumed to 
consist of the alteration of an elliptical orbit resulting from 
satellite launch to a highly circular orbit. This might be re- 
quired in order to eliminate guidance or propulsion errors re- 
sulting from the launching operation. The propulsion system 
is required to convert the excess kinetic energy of the elliptical 
launching orbit to potential energy in order to restore the 
balance required for a circular orbit. 

A comparison of radial and circularizing thrust orientations 
indicate that eccentricity could be reduced more effectively by 
the use of the circularizing thrust pattern. Fig. 4 contains 
the.results of several mission calculations with the circulariz- 
ing thrust pattern. The initial orbit in each case had an ec- 
centricity of 20 per cent and a semi-latus rectum of 25,000 
miles. 

For a (10)~? thrust to weight ratio, thrust initiation at 
apogee appears to result in the minimum time required to 
circularize the orbit, 0.875 hr. It also results in an increase in 
orbit semi-latus rectum. For a (10)~* thrust to weight ratio, 
thrust initiation at perigee resulted in an increase in the semi- 
latus rectum and the smaller propulsion time requirements. 
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Altitude Correction Mission 


‘he altitude correction mission shall be assumed to consist 
of the alteration of an initial low altitude-near circular orbit 
to a higher altitude orbit of the same eccentricity. This 
would conceivably be necessary in order to eliminate launch 
errors or to augment booster capabilities. The propulsion 
system is required to increase the potential energy of the orbit 
and to remove any eccentricity developed during the process. 

lig. 5 contains the results of a comparison between the con- 
stant eccentricity and the transverse thrust orientation pat- 
terns for a thrust to weight ratio of (10)~*. Note that the use 
of the transverse pattern must be followed by a period of cir- 
cularizing thrust orientation. Initial semi-latus rectum was 
4301 miles and terminal value 26,000 miles. 

The desired altitude is achieved at constant eccentricity in 
18.25 hr. The alternate program involves the use of trans- 
verse thrust for 9.5 hr, a coasting period of 5.8 hr and the use 
of circularizing thrust for 1.3 hr. This results in a total pro- 
pulsion period of 10.8 hr. 


Summary 


The precision adjustment of space vehicle orbit characteris- 
tics can be accomplished most effectively by the use of thrust 
orientation patterns peculiar to the requirements of the 
mission. Propulsion time requirements for converting a low 
altitude orbit to a parabolic escape orbit can be minimized by 
the use of tangential thrust orientation. Eccentric orbits can 
be converted most effectively into circular orbits by the use of 
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the circularizing thrust pattern. Altitude corrections can be 
accomplished most effectively by the use of transverse thrust 
orientation, a coasting period and a period of circularizing 
thrust orientation. 
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Minimization of Characteristic Velocity 
for Two-Impulse Orbital Transfer 


H. MUNICK,! R. McGILL? and G. E. TAYLOR® 


Grumman Aircraft Engineering Corp., Bethpage, N. Y. 


Analytic solution is given to the problem treated by 
Lawden, that of optimum transfer from a point on an 
elliptical orbit to a circular orbit. More generally, the 
absolute minimum characteristic velocity path is given for 
transfer between two arbitrary terminals. 


N 1925 Hohmann (1) first established the minimum char- 
acteristic velocity path between circular orbits. In recent 
years Lawden (2) gave a suggestive treatment to the more 
general problem of optimum transfer from a point on an 
elliptical orbit to a circular orbit. The following problem 
was first formulated by Vargo (8). 

The problem treated here is that of finding the minimum 
characteristic velocity path between two co-planar terminals, 
and from a terminal to a “point.”’ For the purposes of this 
paper a terminal is defined as a locus specified by a radial 
distance and a velocity vector. A “point” is specified by a 
velocity vector and a radius vector. The class of orbits con- 
sidered is that for which the apogee of the ellipse correspond- 
ing to the initial terminal is less than the perigee of the ellipse 
corresponding to the final terminal. This is tantamount to 
saying that a minimum of two impulses is required to accom- 
plish the transfer. 

The formulation here essentially follows that of Lawden 


(2). 
Analytical Treatment 


Consider the problem of optimum transfer of a space vehicle 
between two terminals in an inverse-square field. Let (u, 
%) be the velocity components of the first terminal at a dis- 
tance a from the attractive center. Let (ur, vr) be the veloc- 
ity components of the second terminal at a distance 8 from 
the attractive center. At the first terminal an impulse is 
applied resulting in new velocity components (uw, v:). The 
space vehicle then goes into a transfer orbit arriving at the 
second terminal with velocity (ue, v2). Upon arrival at the 
second terminal a second impulse is applied adjusting the ar- 
rival velocity to the desired terminal velocity (ur, vr). For 
such a maneuver the characteristic velocity is given by 


A* = V( ur — wo)? + (01 — 00)? + V (us — up)? + (02 — oF)? [1] 
From conservation of angular momentum 
uz = (@/B)w [2] 


and from conservation of total energy 


v2? = (: + + Qu = [3] 


Introduce dimensionless parameters z, y by 


Ui hal [4] 


Let the distance ratio r be defined by 
r= a/B [5] 


After nondimensionalizing,® the physical problem of mini- 
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mizing characteristic velocity becomes equivalent to the 
mathematical problem of minimizing 


Mx, y) = V(x — a0)? + (y — wo)? + 
V (re — zp)? + — + y? + Ar — 1) — [6] 


subject to 
= (1 O<r<l 


The domain of the independent variables (zx, y) is then a 
region bounded by the ellipse given by Equation [7] and 
circle of arbitrarily large radius (Fig. 1). 

Now consider the region bounded by the ellipse yz = yr and 
a circle of arbitrarily large radius (Fig. 2). For this region, 
from Equation [6] we have 


Xz, y) = V (a — 2)? + (y — yo)? + V (rx — xr)? = X(a,y) (8) 


At an interior point in the region bounded by the ellipse y2 = 
ye and a circle of arbitrarily large radius (Fig. 2), A; cannot 
have a relative minimum. This is so since the necessary con- 
dition that the partial derivatives of \; with respect to 2 and // 
vanish simultaneously requires that r? = 1, which is a contra- 
diction. 

Special consideration is given to the coordinate (2%, yo) anc 
points along the line x = zr/r. The point (2, yo) must lic 
inside the boundary curve yz = 0, otherwise the transfer 
could be made using one impulse. The points along the line 
x = £r/r will be treated subsequently. 

Excluding the point (xo, yo) and points on the line x = 
xr/r from the region under consideration insures that both: 
partial derivatives of \, exist at all points in the region. The 
function , (x, y) is continuous in and on the region. Since 
this region is closed, the function ), (x, y) will assume its least 
value somewhere in this region. However, since the function 
has no relative minima inside the region, it must assume this 
minimum value on one of the boundary curves. On the 
arbitrarily large outer boundary curve, we observe from 
Equation [8] that the function becomes arbitrarily large. 
The conclusion is that \; (x, y) assumes its least value on the 
curve y2 = Yr. 

For the region bounded by the curve yz = yr and a circle of 
arbitrarily large radius, the smallest value that the function 
A(z, y) can hope to attain would be at the point where A, (2, 
y) takes its least value. Since this point occurs on the curve 
Yo = yr, A(2, y) also takes its least value on this curve. 

Now returning to the true domain given in Fig. 1, we have 
established that the minimum of A(z, y) must lie in or on the 
annular region bounded by the curves y2 = 0 and y2 = yr 
(Fig. 2). This drastically reduces the number of possible 
candidates for the solution. 

Investigating the partial derivative of \(x, y) with respect 
to x we establish that this derivative is strictly positive along 
the line x = zr/r. Therefore, since this derivative can never 
vanish along this line, points on the line are not valid candi- 
dates for a relative minimum. 

Examining the partial derivatives of A(z, y) on the curve 
Yo = Yr we observe that for both derivatives to vanish simul- 
taneously r must be equal to + 1. Since this is a contradic- 
tion we conclude that no relative minima exist on the curve 
Y2 = Yr. 

Approaching the boundary curve y, = 0 along rays parallel 
to the x axis we observe that in the immediate neighborhood 
of the boundary curve, OA/Oz ~ — ~. Therefore approach- 
ing the boundary curve yz = 0 along rays parallel to the x 
axis forces A(x, y) to increase. We conclude that on the boun- 
dary curve y2 = 0, A(z, y) cannot assume its least value. 

From these arguments it is finally concluded that the func- 
tion A(z, y) must assume a relative minimum at an interior 
point of the region bounded by the curves y2 = 0 and yz = 


Yr. 
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Solution 


We have been studying the mathematical problem of mini- 
mizing a nonlinear function of two variables subject to an 
inequality constraint. The foregoing analysis has proved the 
existence of an absolute minimum for the function at an inter- 
ior point of a closed region. At this stage we are permitted to 
apply the necessary conditions for a solution, namely that 
both partial derivatives of the function A(z, y) vanish simul- 
taneously. 

Setting both partial derivatives to zero, we deduce that a 
solution must satisfy 


rip? — + + (yo! — 1) V2/ + Xx 
[zeye’ — (we 1) + =0 19) 
where yo’ = Yo/yr. Factoring yields 
fre — [ + (ye — 1)V 2/1 + X 
{re — — — 1) V2/(1 + = 0 [10] 
From the first factor in Equation [10] and the expressions for 


the partial derivatives OA/Oz, OA/Oy we conclude 


YP [11] 


= Yo 
= [12] 
m+rV2/(1+7r) 


Equation [12] gives y as a linear function of x, and setting this 
linear function of z into Equation [11] yields a quadratic 
equation in x whose coefficients depend on the given data. 
We then solve analytically and determine two exact values of 

Using the second factor in Equation [10] and proceeding as 
before we obtain 


YF 
= [13] 
ro— pp — V2/(1 +1) 


yY — Yo Yo 
14} 
2% gy —rV2/(1 +1) 


Similarly, the pair of equations [13 and 14] can be solved 
analytically yielding two more exact values of x. Using 
Equations [12 and 14], we determine the corresponding y 
values. 

We have now determined analytically four solutions (2, y) 
to the minimum problem. For these four values of (x, y) 
there correspond four values of \. The smallest of these four 
values of A is the absolute minimum for \. This then gives 
the complete solution to the terminal-to-terminal problem 
first formulated by Vargo (8). 
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Nomenclature 


u normal component of velocity 

v radial component of velocity 

x = nondimensional normal component of velocity just after 
first impulse 


y = nondimensional radial component of velocity just after 
first impulse 

a = initial radial distance 

8 = final radial distance 

r = distance ratio, a/B 
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Fig. 1 True domain 
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Fig. 2 Annular region 


X = nondimensional characteristic velocity 

* = characteristic velocity of transfer 

mw = gravitational constant 

x2 = nondimensional normal component of arrival velocity 

ye = nondimensional radial component of arrival velocity 
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Pierce Gun Design for an Accelerate- 
Decelerate Ionic Thrust Device 


MANFRED J. RAETHER! and ROBERT N. SEITZ? 
Army Ballistic Missile Agency, Redstone Arsenal, Ala. 


ITHIN the last few years, considerable interest has 

developed in the possibilities of ionic propulsion systems 
for low thrust space propulsion applications. Systems and 
missions analyses have indicated that specific impulse values 
ranging from 5000 to 10,000 sec will be desirable for ionic pro- 
pulsion systems in the near future. For heavy ions, this cor- 
responds to acceleration voltages ranging from about 1.6 to 
6.4 kv. The current densities that can be obtained at these 
voltages for technically feasible acceleration structures are 
rather low (of the order of 1 to 3 ma/em?). Consequently, 
large emitting areas or a large number of ion sources are re- 
quired for high currents. This, in turn, results in a consider- 
able increase in weight, volume and complexity. 

These difficulties may be avoided to a certain extent by ex- 
tracting the ions from the ion source with a high voltage and 
decelerating them to the desired lower voltage. If there is no 
interception of current by the acceleration electrode, no net 
power is consumed during the acceleration-deceleration proc- 
ess. 


Conditions for Effective Deceleration 


We shall now specify the conditions, under which an ac- 
celerate-decelerate technique is feasible. We consider a 
space charge flow in a parallel plane geometry (Fig. 1). 

It will be assumed that the ions are emitted from the anode 
with zero initial velocity. In region (1) they are accelerated 
to the potential V.. They now enter into region (2) with an 
initial energy eV., where they are decelerated to the desired 
potential V,. 

The problem is now reduced to the discussion of space 
charge flow with finite initial velocities in a planar diode. 
This case has been treated by Fay, Samuel and Shockley (1).* 
Their results will be briefly discussed in the following para- 


graphs. 
We introduce the following dimensionless quantities 
= x/s = V/V, 
where 
4 Qe 1/2 
=) 


coordinate in the direction of flow measured from the 

acceleration plane 
= distance between the anode and acceleration planes 

s = distance between the acceleration and deceleration 
planes 

potential 

current density 

= absolute magnitude of the electronic charge 

= jonic mass 

= permittivity of free space = 8.85 X 10-!* coulombs per 
v-m 


z 


Mks units are used throughout this paper. 

If we inject a space charge limited current into the de- 
celeration region, 7 becomes s?/d?. 

The quantity Jo is used as a convenient reference current 
density in the dimensionless ratio 1 = I/Io. 

In Fig. 2, the possible types of space charge flow are shown 
in an (7, @,) diagram. @,, the deceleration ratio, is the de- 
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Fig. 1 Space charge flow in an accelerate-decelerate parallel 
plane geometry 
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Fig. 2. Possible types of space charge flow in the deceleration 
region 


Fig. 3 Diagram of operating conditions for unimpeded space 
charge flow in the deceleration region (assuming space charge 
limited injection current) 
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celeration potential measured in units of the acceleration 
potential. 

In region B of this diagram, the potential has a minimum, 
and dmin = 0; i.e., a virtual anode is formed. Only part of the 
injected current is transmitted. 

In region C, the potential also has a minimum, but ¢min 
> 0. All the injected current is transmitted. In region D, 
no minimum occurs; all current is transmitted. 

In the overlapping regions B, C and B, C, D, all the re- 
spective types of space charge flow are possible, the actual 
flow pattern depending upon previous variations in the de- 
celerating voltage. The shaded region does not permit stable 
operation. 

If it is necessary that all the injected current be trans- 
mitted, we are restricted to the region between the curve c 
and the ¢, axis. It is noted that, theoretically, the ions may 
be decelerated to zero velocity without loss in current. There 
may, however, be a practical limitation upon the deceleration 
ratio. It may be required that the field strength in the de- 
celeration gap must not exceed the field strength in the ac- 
celeration gap. This would mean that 


Veg 
8 d 


or 
s/d > 1 = ¢, 


Fig. 3 shows the region of interest in an (s/d, ¢,) diagram for a 
space charge limited injection current. The choice of a point 
inside the shaded area specifies s/d and the deceleration ratio. 
It is seen, that based upon the foregoing assumption, the 
maximum practical deceleration ratio would be about 1:20. 


Application of the Pierce Gun Design Principles 


The rectilinear space charge flow discussed previously may 
be realized by applying the Pierce gun principle to both the 
acceleration and deceleration region (2) for a rectangular-slit 
gun geometry. [See Reference (2).] Fringe effects at the 
ends of the slit are neglected in this treatment. 

We shall discuss the application of the Pierce gun design to 
the deceleration region for a special set of the parameters s/d 
and ¢,. We choose s/d = 1 and ¢, = 0.1; this point lies 
within the D region. The potential distribution for the D 
region is given by (1) 


where (d/s)'? and C(d/s)!? are the respective argu- 
ments of the function y. wy is defined by 
W(x) = (1 + — 22) 


C is a constant of integration which must be determined from 
the boundary condition that for § = 1, ¢ = ¢, = 0.1. For 
our special example we find 

C = —0.2 ¥(C) = 1.25 


To obtain the electrode shapes, we follow the method out- 
lined by Pierce and replace & by the complex quantity € + jn, 
where j = (—1)!*. We find 


erin = 


where ¢ is now a function of § + jn. can be separated into 
its real and imaginary parts 
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Fig. 4 Electrode configuration for accelerate-decelerate Pierce 
gun (10:1 deceleration ratio) 


This yields 


d\ 
d\ 1/2 d\1/2 
ye + jJ)2+C (‘) ce + 2 (‘) c| 


In our special case, we have 


E+ jn = 125 + + 04] 


This complex equation is equivalent to two real equations 
from which J must be eliminated. The curves R = R(é, 7) = 
constant are then the equipotential lines. Any of these lines 
may be replaced by an electrode. The foregoing equation, 
however, does not allow a separation into its real and imagi- 
nary part in a convenient way. Therefore, it is better to re- 
tain J as a parameter and to present the curves in the para- 
metric representation 


E=ER,J) 1 = (PR, J) 


The inner and outer electrodes are given by 


E=&1,J) = 7(1,/) 
f= %0.1,J) 9 = (0.1, J) 


respectively. 

The equipotential lines have been calculated numerically, 
and the results are shown in Fig. 4, together with the Pierce 
gun design for the acceleration region. 

Between the acceleration region and the deceleration re- 
gion, a field-free gap remains. The length of this gap is com- 
pletely arbitrary. However, the length of this gap, as well 
as the choice of beam width aperture, is contingent upon three 
major types of electrostatic beam distortion. First, distortion 
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arises because of the finite extent of the electrodes, which 
ideally should be infinite. This effect may be minimized 
by making the field-free gap as long as possible. Second, 
distortion arises because of the lateral divergence of the ion 
beam resulting from the mutual repulsion of the ions. This 
beam spreading effect is dependent upon beam width and 
depends also upon the length of the field-free gap. This type 
of distortion may be minimized by making the gap as short 
as possible. A third type of beam distortion is caused by the 
divergent-convergent lens effects resulting from the fringing 
fields at the entrance and exit apertures of the field-free space. 
It may be minimized by reducing the beam width. Thus, in 
choosing the length of the field-free gap, an optimum balance 
exists between the first and second types of distortion. In 
any event, the beam width should be as small as possible. 


Conclusions 

Since a number of simplifying assumptions have been made 
in deriving these results, this treatment cannot be regarded 
as more than a first approximation to a final lens design. 

It should be borne in mind that maximum advantage may 
be realized from an accelerate-decelerate system only if prac- 
tical structural limitations rather than voltage breakdown 
limitations determine the acceleration gap d. The available 
current densities are also limited by the current densities 
which may be obtained from the ion source. 
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Magnetohydrodynamic Cavities 
BERNARD STEGINSKY' 


Battelle Memorial Institute, Columbus, Ohio 


Some features of magnetohydrodynamic cavity forma- 
tion are discussed in qualitative terms. Examination 
of the dependence of the final steady-state configuration 
on the initial conditions in two cases reveals incompat- 
ibilities which may arise in the physical interpretation of 
the problem if a proper ‘‘model’’ is not specified. 


EVERAL theoretical papers have recently appeared deal- 

ing with the existence of cavities in magnetohydrody- 
namics (1—4).2 These cavities arise in connection with the 
flow of a perfectly conducting fluid past a magnetic field source 
such as, for example, a line current or a magnetic dipole. A 
perfectly conducting fluid cannot penetrate a region where 
the magnetic field is too strong, nor will such a fluid “sup- 
port”’ diffusive transport of the field. Consequently, it is 
possible to demonstrate that the considered flow will result 
in the formation of a cavity region surrounding the source, 
around which the fluid flows as it would around a solid body, 
and in which the entire magnetic field strength is confined. 
The shape of the cavity-fluid boundary may be determined 
from the given distribution of current within the cavity, the 
prescribed type of external flow, and the condition of pressure 
equality across the interface. 

In aerodynamic applications, we are usually interested in 
obtaining the steady-state solutions corresponding to a fina] 
equilibrium configuration. That is, we wish to determine 
the equilibrium shape of the cavity boundary subject to a 
given flow, without regard to the nature of the transient for- 
mation process from some initial state. All too often, there- 
fore, the role of the initial conditions in determining the ulti- 
mate shape acquired by the cavity is overlooked. This over- 
sight has already resulted in the appearance of a paradox in 
a recent treatment of the problem. The purpose of the 
present note is to clarify this particular aspect of the problem 
in qualitative terms. 
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In considering initial states of the fluid, we shall distinguish 
between the two limiting cases of zero conductivity and in- 
finite conductivity. In both cases, however, the fluid will be 
assumed to be in an initial state of rest. When the conduc- 
tivity of the fluid is infinite, a magnetic field source will pro- 
duce a cavity region surrounding the source (a circular shape 
in the case of a line current), the interior of which contains 
the entire magnetic field but no fluid. At the cavity surface, 
the hydrostatic pressure outside balances the magnetic pres- 
sure just inside the surface. Cole and Huth (2) have studied 
this problem for the cases where the source is represented by 
a, a line current and 6b, a magnetic dipole (the limit of two 
approaching line currents). They have shown that, when 
an incompressible flow field is superposed on a, the original 
shape is distorted into an ellipse with major axis transverse 
to the direction of flow. The elliptical shape results from a 
compression of the cavity parallel to the flow direction, be- 
cause of the increase in stagnation pressure and deficit in 
hydrostatic pressure at the points of maximum velocity. 
Although only the static solution for b was obtained, Cole 
and Huth indicate that the original shape for this case would 
not be much distorted by a flow. 

The second case, of a nonconducting fluid, is fundamentally 
different from the first one in that there will be no cavity 
region to begin with, and the magnetic field will penetrate the 
whole fluid. We may then inquire as to the final configura- 
tion of the system if, suddenly, the fluid were endowed with 
infinite conductivity and a given motion. This problem, con- 
siderably more difficult than the previous one, has been 
discussed qualitatively in (1 and 6). 

The following idealized example will be used to illustrate 
the interaction between the fluid and the magnetic field lead- 
ing to the formation of a cavity in the flow. Suppose that, 
in a nonconducting fluid which is at rest, there is placed a 
body having in its interior an arrangement of current-carry- 
ing coils. The magnetic lines of force pass continuously 
through the body surface and penetrate the surrounding fluid 
to great distances. 

Now let the fluid be suddenly endowed with infinite con- 
ductivity and a general motion past the stationary obstacle. 
With the start of the motion, in which all parts of the fluid 
participate, the magnetic lines of force will be swept down- 
stream by the fluid. This follows from the theorem of 
magnetohydrodynamics which states that magnetic lines of 
force move about in a perfectly conducting fluid as thoug) 
rigidly attached to the fluid particles. This sweeping motion 
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results, however, in compression of the lines of force with a 
corresponding increase in magnetic pressure. The magnetic 
pressure tends to deflect the fluid away from a region of strong 
magnetic field concentration. Alternatively, the decay time 
of the component of fluid velocity transverse to magnetic 
lines of force varies inversely with the fluid conductivity. 
Therefore, motion across the lines of force will be annihilated 
almost instantaneously. 

The system will very rapidly equilibrate to a state wherein 
the following conditions are satisfied : 

1 A cavity region forms around and behind the body, ex- 
tending far down stream, through which the fluid cannot 
penctrate; the fluid velocity is tangent to the surface of the 
cavity. 

2 The magnetic field is confined to this cavity region, 
since the field cannot penetrate a perfectly conducting fluid; 
the magnetic lines of force are tangent to the surface of the 
cavity, and a current sheet of infinite density forms on the 
surface. 

3 The external hydrostatic pressure balances the mag- 
netic pressure across the cavity-fluid interface. 

‘he above conditions, together with the given distribution 
of current in the body and a prescribed type of flow, are 
suflicient to determine the shape of the cavity. V. N. 
Zhigulev (3), making clever use of complex variables, recently 
solved the two-dimensional problem of this type for the flow 
past a line current. The flow was chosen to correspond to a 
Newtonian pressure distribution at the cavity surface. 

‘The final equilibrium configuration could conceivably be 
achieved by the passage of a high speed shock wave through 
the initial system, as suggested by Burgers (1). The gas 
immediately behind the initial shock wave moves at a hyper- 
sonic speed which is slightly less than the propagation ve- 
locity of the wave. Therefore, a stationary bow shock wave 
forms ahead of the body. The hot, ionized gas generated 
by the initial shock wave possesses a relatively high conduc- 
tivity and can therefore interact with the magnetic field 
emanating from the body. 

Because the electrical conductivity of the fluid actually is 
finite in practice, the entire cavity region would be occupied 
by fluid in a state of rest or near-rest since leakage of the fluid 
across the magnetic force lines is possible. This region, where 
the essential strength of the magnetic field is concentrated, is 
the so-called magnetic boundary layer. 

Cavity shapes in perfectly conducting fluids may be ob- 
tained to a good approximation from solutions of magnetic 
boundary layer problems for very large, but finite, magnetic 
Reynolds numbers. Ludloff (5) has recently shown that the 
flow of a fluid of large magnetic Reynolds number past a 
line current results in the formation of a thin magnetic 
boundary layer which is greatly elongated in the flow direc- 
tion. This result is markedly different from that of Cole and 
Huth for the case of a perfectly conducting fluid, in which an 
originally circular cavity experiences an elongation in the 
transverse direction and a contraction in the longitudinal 
direction. 

This seemingly anomalous behavior may be resolved by 
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observing that the superposition of a flow on the two initial 
states produces a fundamentally different effect in each case, 
leading to dissimilar end results. In the first case, the con- 
ductivity of the fluid was assumed to be infinite prior to the 
onset of motion. Thus, the fluid cannot interact directly 
with the magnetic field, which is zero everywhere in the fluid, 
but distorts the original shape of the cavity by virtue of the 
redistribution of pressure on the boundary. In the second 
case, the conductivity of the fluid was considered to be zero 
in the rest state and only assumed an infinite value simultane- 
ously with or after the commencement of flow. Therefore, 
in this case, the fluid can interact directly with the magnetic 
field and, in doing so, sweeps the magnetic force lines down- 
stream. 

In a recent paper, Sakurai (4) studied the problem of the 
“steady, two-dimensional hypersonic flow of an ideal gas with 
infinite electrical conductivity around a two-dimensional 
magnetic dipole, the axis of which is perpendicular to the 
direction of the uniform flow.’”’ The cavity region in the 
steady-state (or “frozen region,’ as Sakurai called it) was 
assumed to have a circular shape, and the magnetic field 
within the cavity was represented by the solution for the 
static magnetic dipole. 

From a mathematical standpoint, the problem becomes over- 
determined if, in addition to the condition of pressure equality 
at the interface, a given pressure distribution in the external 
flow and a given magnetic source distribution within the 
cavity, we also attempt to specify the shape of the cavity. 
Consequently, it must be assumed that the cavity shape is 
unknown and is to be determined from the conditions of the 
problem. 

In the light of the previous discussion, it would seem 
plausible to conclude that the steady-state cavity region could 
not assume as compact a shape as that of a circle unless the 
conductivity of the fluid were infinite initially, i.e., before the 
commencement of flow. In the opposite limiting case of a 
zero conductivity initially, we would expect the formation of 
an extended cavity region in the downstream direction. Al- 
though neither of the two cases examined is an exact repre- 
sentation of hypersonic flight conditions, the latter furnishes 
the more accurate “model” for the problem. It would then 
appear that Sakurai’s specification of a circular cavity, apart 
from its inherent mathematical difficulty, is also incompatible 
with the conditions of hypersonic flight. 
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Analytic Formulation for Radiating 
Fins With Mutual Irradiation’ 


E. R. G. ECKERT,? T. F. IRVINE Jr.* and 
E. M. SPARROW? 


Heat Transfer Laboratory, University of Minnesota, 
Minneapolis, Minn. 


HE ABSENCE of convection in the operating environ- 

ment of a space vehicle has directed attention to radiation 
as the primary or exclusive transferring mechanism for dis- 
posing of waste heat. Radiative collectors are also being 
proposed for capturing solar energy for power generation 
purposes. The use of finned surfaces, already a well-known 
technique in convection heat exchange equipment, is now 
being considered in the design of space vehicle radiators. The 
fins would not only augment the heat transfer characteristics 
of the radiator, but would also provide protection for the base 
surfaces against bombardment by space particles. 

In studying the characteristics of a radiating-conducting 
fin, the first step which has been taken (1—4)‘ is to suppose that 
there is no radiant interaction with adjacent fins or with the 
base surface. For such an isolated fin, the conservation of 
energy principle can be expressed as a differential equation, 
and results have been obtained for both the plane fin and the 
circumferential fin of constant thickness. In almost all real 
situations, there will be a radiant interchange between ad- 
jacent fins and with the base surface as well as with the ex- 
ternal environment. Part of the radiation leaving a given fin 
will be intercepted by a neighboring fin or by the base surface, 
and, if the surfaces are nonblack, there will be reflections and 
re-reflections. The mathematical description of this complex 
interchange process, which includes both energy transfers 
from physically separate surfaces and heat conduction, must 
be made as an integro-differential equation. The formulation 
of such problems may be unfamiliar to many heat transfer in- 
vestigators, and it is the purpose of this report to show how 
such an analysis is carried out. Prior consideration of the 
radiant interference between fin and base surfaces has been 
made in (3 and 5), but in an incomplete and restrictive 
manner. 

Before proceeding to the mathematical formulation, it is 
worth while to mention some general characteristics of radiat- 
ing fins. Depending upon the radiation and conduction prop- 
erties of the material and on geometrical parameters and 
orientation, the use of radiating fins may result either in an 
increase or in a decrease in heat transfer. For example, con- 
sider a plane wall which is of infinite length in the direction 
normal to the page [Fig. 1(a)]. For finite fin conductivity and 
black surfaces, the presence of fins definitely decreases the 
total heat transferred, relative to the unfinned situation. On 
the other hand, if the emissivity is less than 1, then it can be 
said that fins of infinite conductivity will definitely increase 
the total heat transfer, and that increases are also likely for a 
range of finite conductivities. As a second illustration, con- 
sideration is directed to the long, longitudinally finned cylin- 
der shown in Fig. 1(b). Here, in contrast to the previous ex- 
ample, finning may increase the heat transfer even for finite 
fin conductivity and black surfaces. For nonblack surfaces, 
the remarks made in connection with Fig. 1(a) continue to 
apply. Thus, from this, there is indicated the need for careful 
analytical study prior to the use of radiative finned surfaces. 
It is evident that there will be a large number of parameters 
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(a) Longitudinal fins on 
a plane surface 


(b) Longitudinal fins on 
a cylinder 


Fig. 1 Some fin configurations 


which will determine the heat transfer performance of radi.:t- 
ing fins, more than are usually found in the convective ca-e. 
Consequently, analysis of finned radiating surfaces may in 
many cases have to be carried out with a specific application 
in mind. 

With this as background, we now proceed to describe the 
mathematical formulation for ensembles of radiating-condu  t- 
ing fins. Attention will first be focused on a fairly general ‘in 
configuration; then, the governing equations will be given tor 
some specific cases for which numerical calculations are now 
under way. 


General Mathematical Formulation 


The analytical approach to be presented here applies to a 
variety of geometrical configurations. But, for the sake of 
concreteness, the derivation will be carried out with referen:e 
to an easily visualized physical situation. Consideration is 
given to the conduction-radiation transfer process for a long 
circular cylinder which is finned either longitudinally or cir- 
cumferentially as shown in Fig. 2, in (a) and (b) respectively. 
The longitudinal fins are assumed to be arranged symmetri- 
cally around the cylinder, whereas the circumferential fins are 
equally spaced along the length of the cylinder. Only a 
typical pair of neighboring fins is shown in the figures. The 
thermal conditions within the cylinder and in the external 
environment are uniform. The radiation coming from the 
environment is characterized by an equivalent black body 
radiation e (per unit area) passing uniformly through the area 
S, connecting the tips of the fins (dashed line in Fig. 2). The 
surfaces of the fins and of the cylinder are taken as diffuse 
emitters and reflectors. The surface of the cylinder, i.e., the 
fin base surface, has a uniform temperature 7>. Initially, it 
will be assumed that the surfaces are gray (a = € = 1 — p), 
and also, that the fin surfaces are convex or plane, so that no 
surface element of a fin ‘“‘sees”’ any other part of the same fin. 
The modifications associated with lifting these assumptions 
will be discussed later. 

The starting point of our study is the basic principle of 
energy conservation. For steady-state conditions, the appli- 
cation of this law to the cross-hatched element of Fig. 2 yields 


(Qnet Jeond (Qnet )rad =0 [1a] 


The conduction term is evaluated in the usual fashion using 
Fourier’s law, giving 


To characterize the net radiation, it is of advantage to intro- 
duce two new quantities. [See (6), pp. 416-418.] First, 
there is B, which represents the combined radiant flux 
(emitted plus reflected) which leaves a position on the surface 
per unit time and area. The second quantity is H, which is 
the incident radiant flux at a position per unit time and ares. 
With these definitions, the net radiant flux leaving the shaded 
element is 


H(zx)| [2] 
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(a) Longitudinal fins 


(b) Circumferential fins 


Fig. 2 Diagrams used in derivation of governing heat tranfer equations 


From the geometry of Fig. 2, dS, = l(x)dz/cos a, where I(x) 
is the fin perimeter. With this, the conservation of energy 
principle [1b] becomes 


d Ea (B(x) — H(2)) [3] 


dx dx cos @ 


As it stands, Equation [3] contains three unknowns: T(z), 
B(x), and H(z). So, it is necessary to invoke other relation- 
ships between these quantities. First, we note that the com- 
bined radiant flux B leaving a surface is composed of two 
parts: Emitted energy and reflected energy. The emission is 
given by the usual expression 


eoT (x) 
whereas the reflection is represented in terms of the incident 
energy and the reflectivity as 
pH(x) = (1 — €)H(z) 

With’ this, the definition of the combined flux B can be 
evaluated as 

B(x) = eo T(x) + (1 — €)H(z) [4] 
The flux H which impinges on z is related to the energy com- 
ing from the external environment and to the flux leaving the 
other surfaces, i.e., the opposite fin and the base surface. 
Considering first the contribution of the opposite fin, we note 


that the radiant energy leaving a typical area dS, in all 
directions is 


B(y)dS, [5a] 
Of this, an amount 
B(y)dS,dF [5b] 


arrives at x, where dF,_, is the angle factor [see (6), pp. 396- 
398] under which dS, is seen from y. Using the reciprocity 
relation dS,dF,_. = dS.dF,_,, expression [5b] becomes 


B(y)dS.dF ,_, [5e ] 
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So, the energy impinging per unit area at dS, due to radiation 
leaving dS, is 


B(y)dF [5d] 


But, dS, receives energy from all positions y, and hence, the 
total contribution of the opposite fin is found by integrating 
expression [5d] giving 


ffi, (6 


In the same way, the contributions of the base surface and of 
the external environment to the incidence at dS, can be 
formulated. With these, H(x) can be written as 


H(x) = B(y)dFr-y + eF [7] 


where ii is to be remembered that dF,—, is zero whenever x 
and y cannot see each other. Because of the symmetry, it is 
apparent that the distribution of the combined flux B over 
the right-hand fin will be the same as the distribution over the 
left-hand fin. So B(x) and B(y) are the same function; only 
the independent variables have been interchanged. 

Summing up thus far, we find that it has been possible to 
augment the energy conservation Equation [3] with Equations 
[4 and 7] which essentially express conservation of radiant 
flux. But, in this process, the previous unknowns 7x), B(x) 
and H(z) have been joined by the radiant flux B(z) of the 
base surface. To make the problem determinate, additional 
flux conservation equations can be written for the base sur- 
face in a manner analogous to Equations [4 and 7]. The en- 
tire set of governing equations is [3, 4 and 7] and 


B(z) = eoT,! + (1 — €)H(z) [8] 
H(z) =2 f, + 19] 


To complete the statement of the problem, the following 
boundary conditions are prescribed: 


When zx = 0 


dSy | d S x ae, 
A(x) 
dx 
e | 
\ ds dS 4 
y 
\ 
A(x) | | 
\ 
pe 
eS 
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When z = L 
dT(x) B(L) H(L) 
dz sin @ 


or 
dT(x)/dz = 0 


The second condition at z = L is usually taken as an approxi- 
mation of the first one. 

Once the geometry of the system has been specified, the 
angle factors F may be evaluated. Then, the set of five 
nonlinear integro-differential equations may be solved to ob- 
tain all heat transfer quantities of interest. For example, 
when B(z) and H(z) have been determined, the local heat loss 
may be found from Equation [2]. The solution of the govern- 
ing equations will almost certainly have to be carried out 
numerically. 

To extend the equations to fins with concave surface, it is 
only necessary to add to the right side of Equation [7] a term 
which accounts for the irradiation of the element dS, by other 
parts of the surface of the same fin. For fins whose surfaces 
have absorptivities and emissivities which depend on wave 
length, Equations [4, 7, 8 and 9] are written for monochro- 
matic radiation. Then, in Equation [3], B and H have to be 
introduced after integration over all wave lengths. The actual 
mechanics of carrying out such a calculation appears ex- 
tremely involved. 

For the special case of black surfaces (€ = 1), the first three 
governing equations can be solved without recourse to the last 
two. After rearrangement, Equations [3, 4 and 7] yield 


d aT(z) 
Ee 


U(x) 
COs @ 


This same derivation is also valid if pins, rather than fins, 
are used as the extended surfaces. It is only necessary to re- 
interpret the integral over S, in Equations [7 and 12] as in- 
cluding the surface area of all of the pins which may be seen 
from the typical pin S,. A similar interpretation is to be 
made of the integral over S, in Equation [9]. 


Some Special Cases 

The configuration of Fig. 2 has been useful for demonstrat- 
ing the method of formulating the radiative fin problem, but 
it contains too many parameters to be an appropriate case for 
numerical study. Consequently, consideration is being given 
to two somewhat simpler configurations, Fig. 3(a) and (b), 
which still have practical utility and whose solutions should 
provide insight into the behavior of radiating fin ensembles. 
Fig. 3(a) represents a typical pair of longitudinal fins belong- 
ing to a group which are symmetrically arranged around the 
circumference of a long cylinder. As a simplification, the 
curved base surface has been neglected, but this should have 
small influence as long as the fin length L is markedly larger 
than the cylinder diameter. In Fig. 3(b), there are depicted 
two members of a set of circumferential fins which are ar- 
ranged along the length of a cylinder. Once again, the base 
surface has been neglected. 

The governing equations are obtained by specializing Equa- 
tions [3, 4, 7, 8 and 9], which were derived in the previous 
section. Only the final dimensionless results are given here. 
For the radial fins of Fig. 3(a), the derivation yields 


d6(£) 
dé? 


= N,[B*(£) — H*(&)] [13] 


5 For surfaces of infinite extent normal to the page, dF,_, = 
(1/2)d(sin ¢), where ¢ is the angle between the normal at dS, and 
the connecting line from dS, todS,. dF,z, F:_s,, etc. 
are found in the same way. [See Jakob (7), Eqs. 31-49.] 
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&n(1 — cos? 6) 


1 
H*(i) = a) B*(n) (£2 + — cos 6)32 


N, E cos 6 — & 
2 (1 + & — 2€ cos 6)1/2 [14] 
B*(&) = + (1 — €)H*(é) {15] 
where 
B* = B/oT,! 
H* = H/oT,! 
6 = T/T, 
= 2/L 
= y/L 
N, = (L*/kt)oT,? 
N, = e/oT,* 


For the circumferential fin ensemble, Equation [15] remains 
unchanged, and Equations [13 and 14] are replaced by 


1d | , _ — He 
Al N.[B*(é) — H*()] [17] 


ny? + & + 9) 
| 
where y = h/L andr,;* = r;/L <1. The boundary condi- 
tions on both problems are 


= 1 


1 
H*(E) = dn + 


d6/d&(1) = 0 


Inspection of the governing equations shows that the solutions 
will depend on three parameters: The conduction parameter 
N., the environmental irradiation parameter N,, and a geo- 
metric parameter—either the angle 6 or the aspect ratio 
y =h/L. Numerical calculations are already under way, but 
there are no noteworthy results to report at this time. 
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Low Thrust Correction of Orbital 
Orientation 


L. RIDER! 


Convair (Astronautics), San Diego, Calif. 


A method for performing small corrections of orbital 
orientation is described. It is found that the orientation 
parameters of inclination angle and nodal longitude can 
be varied by a low thrust device continuously applied per- 
pendicularly to the orbital plane with sense of application 
reversed twice in every revolution. 


N A RECENT paper (1),? Lorell and Lass have determined 

that the effect of a constant radial perturbing force on the 
orbit of a satellite is to cause the line of apsides to advance 
while maintaining the apogee and perigee distances constant. 
In terms of our notation, the rate of advance of the line of ap- 
sides reduces to 


2rAV 1 — e? radians per revolution (1) 


Derivation of Equations 
Moulton (2) gives the rates of change of 7 and Q due toa 
normal disturbing acceleration W as 


di r cos (v + w) 


= Ww 2 
dt na? VE = 
dQ r sin (vy + w) 
—= Ww [3] 
dt — sini 


in standard notation. 
Conservation of angular momentum provides the relation 


[4] 


Substitution of Equation [4] in [2 and 3] and conversion to ‘ 


the present notation provides 
di 
dy ( 


A(1 — sin (v + w) 


cos (v + w) 


(1 + € cos 


= [6] 
wiere \ is the ratio of the perturbative accleration (assumed dy sini (1 + € cosy)? 
= Equations [5 and 6] can be integrated as 
(1 — é?) J 1-—eé\ 
S di cose ( ) sinw] + 
(1 + 2e?) sin cos w 3€ cos w =)" (7) 
2(1 + € cos v) (1 — €?)!/2 l+e 2 
[sin vin + ( € ) cos | + 
2) gi i si 
(1+ 2¢*)sinvsinw tan| Al (8) 
2(1 + €cosv) (1 — €?)!? 1+e 2 
In terms of eccentric anomaly E, these equations are 
cos sin w 
di =X ((2 — cos E + 2?) sin — FE] — (1 — € cos [9] 
J 
S ((2 — € cos E + 2e?) sin E — (1 — € cos nyt 
[10] 


constant) to the gravitational acceleration at a distance from 
the center of Earth equal to the semimajor axis of the per- 
turbed orbit. 

In the present note, a method employing low thrust is in- 
vestigated for the purpose of changing two other convention- 
ally chosen parameters which describe the orientation of any 
orbit: Inclination and nodal angle. This method in conjunc- 
tion with that found by Lorell and Lass provides a relatively 
simple way to adjust the orientation of any orbit using a low 
thrust device. 

Assumptions for the analysis are a constant perturbing 
acceleration always applied normal to the orbital plane and 
those of a simplified two-body problem. 
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Determination of the Thrust Program 
Maximum orientation parameter changes per revolution 


If Equations [7 and 8] or [9 and 10] are evaluated over a 
complete revolution in an arbitrary orbit and if the thrust 
(perpendicular to the orbital plane) has a programmed sense 
as shown in Fig. 1, we obtain the maximum change in either of 
the orientation parameters as shown in Fig. 2. Generally, 
such a maximum change in one of the parameters will be as- 
sociated with a nonzero change in the other parameter. 

It is of interest to note from Fig. 2 that when w = 7/2, 
Ai = 4A, i.e., the maximum change in inclination per revolution 
is independent of the shape (eccentricity) of the elliptic orbit. 

Similarly, for a = 0, AQ = 4)/sin 7 and is independent of 
orbital eccentricity. 


Single orientation parameter changes 


Circular orbits have the unique property that when the 
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dt r? 
a 


sense of the perturbing acceleration is such that maximum in- 
clination change occurs (see Fig. 3), no nodal change will 
occur. Alternately, when maximum nodal change is pro- 
grammed, no inclination change occurs. 

When w = 0 and thrust is programmed for maximum in- 
clination change in any elliptic orbit, no nodal change will 


occur. 

sabeboenses THRUST APPLIED OUT OF THE PLANE OF THE PAPER 

eee = THRUST APPLIED INTO THE PLANE OF THE PAPER 
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Fig. 1 Sense of normal thrust for maximum orientation changes 
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Fig. 2 Maximum inclination or nodal change per revolution 
648 


When w = 2/2 and thrust is programmed for maximum 
nodal change in any elliptic orbit, no inclination change will 
occur. 

In general, for a single orientation change of particular 
magnitude (or a change of both orientation parameters of 
particular magnitudes), an interrupted thrust program and/ro 
a new location of the thrust change-over points must be de- 
termined. Equations [7 and 8] or [9 and 10] can facilitate 
the calculation of such a program. 


Nomenclature 


i orbital inclination 


Q = longitude of ascending node 

« = perigee angle measured in the plane of the orbit from the 
ascending node to perigee 

Ai = change in inclination per orbital revolution, radians jer 
revolution 

AQ = change in nodal longitude per orbital revolution, radi:ns 
per revolution 

» = true anomaly 

E = eccentric anomaly 

a = semimajor axis 

Jo = u*/a*, acceleration of gravity at a distance equal to a 

W = small constant perturbative acceleration applied perpcn- 
dicularly to the orbital plane 

h = W/ge 

uw? = Earth gravitational constant, GMe 

c = angular momentum per unit mass 
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Some New Satellite Equations 


RAIMOND A. STRUBLE! 
North Carolina State College, Raleigh, N. C. 


ECENTLY the writer has derived completely general 

equations of motion of a satellite relative to the true 
orbital plane, i.e., the plane which at each instant contains 
the satellite and its velocity vector. For the familiar gravi- 
tational potential 


2 4 
= gk — cos? {1] 


r r3 


they reduce (without approximation) to the third order system 
(sec Fig. 1 for notation) 
a JgR*u cos* i sin 7 sin 28 (2] 
+ JgR*u cos i(1 — cos 28) 


d du 
k 78 (i 3) + usec?i = 


gk? JgRtu? 


(1 — 3sin?isin? [3] 


where 


p = angular momentum (constant) about the polar axis 
k = sec i + (2/p?) JgR*u cos* sin? B 
u=1/r 


It should be noted that the argument of the line of nodes 
Q does not appear at all, and hence the radial motion (orbit 
shape) is independent of the precessional motion of the orbital 
plane, a fact well-disguised by recent (approximate) satellite 
theories. On the other hand, the variations of 7(the angle of 
incidence of the orbital plane) are of first importance to the 
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radial motion. An exact orbit (corresponding to [1]) is pre- 
scribed by the system [2 and 3], and once it has been obtained, 
the (exact) precession of the orbital plane is given by a 
quadrature of 


dQ JgR*u cos? i (1 — cos 28) 


dg p? + JgR* u cost (1 — cos 28) [4] 


From [4] we see that the rate of regression of the line of 
nodes is not uniform but varies (in each half revolution of the 
satellite) from zero to approximately twice the familiar mean 
rate. 

Several facets of these equations which concern satellite 
engineering are to be developed in future papers. 


Simple Method for Approximating the 
Characteristics of Low Thrust 
Trajectories 


PAUL D. ARTHUR,'! HANS K. KARRENBERG! and 
HAROLD M. STARK? 


Systems Corporation of America, Los Angeles, Calif. 


OW THRUST trajectories are of interest for interplane- 
tary travel. Unfortunately, general solutions for these 
trajectories do not exist in closed form, and simple approxi- 
mate solutions are needed. The method of trajectory approxi- 
mation presented here is similar to perturbation methods 
standard to celestial mechanics. Only deviations from a ref- 
erence ballistic (two-body) orbit are calculated. The first 
approximation to these deviations is the expected s = 4 X 
at?. However, this expression is not sufficiently accurate for 
most purposes. More accurate approximations are derived 
in this paper. 


Received Feb. 9, 1960. 
1 Member Scientific Staff. Member ARS. 
2 Member Scientific Staff. 


Juty 1960 


Derivation of Equations 


In the following presentation a simplified model of the solar 
system as discussed in (1),3 is assumed. The orbits of the 
planets are assumed to be co-planar concentric circles about 
the sun. During its interplanetary travels the vehicle is as- 
sumed to be solely under the gravitational influence of the 
sun’s central force field. In addition, any thrust vector is 
assumed to lie in the plane of the planets’ orbits. Thus, the 
motion is two-dimensional. 

The vehicle acceleration A due only to the thrust vector T 
is aligned at some angle @ with respect to the radius vector 
from the sun as indicated in Fig. 1. A, and A, are the radial 
and circumferential components of A. The position of the 
vehicle at any time ¢ is defined by a radius vector R and a 


polar angle ¢. 
The equations of motion are 
-R = (*) + A, [1a] 
3 Numbers in parentheses indicate References at end of paper. 
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Fig. 2 Additional parameter definitions 
where 
Ry = 1 astronomical unit 
9 = gravitational acceleration of the sun at Ro 
Changing variables, as in (2, 3 and 4) 
Ro Jo Ro dr V qRo 
Equations [1] thus become 
1 
+e [3a] 
wy) [3b] 
r dt 


where the dots indicate differentiation with respect to r. 

This ballistic (two-body) orbit which osculates with the 
actual trajectory at tr = 0 will be used as a reference orbit. 
The equations of motion of this reference orbit are (subscript 
b denotes ballistic orbit) 


1 


Fo — To (do)? [4a] 
Tb 
1a integrates to 
— = 0 —— rip =H [4b] 


re dr 
where H is the constant angular momentum of the reference 


orbit. 
Equations [3 and 4] may be combined as follows: 
1 Define X * r — ry and substitute into Equation [3]. 
2 Subtract [4a] from [8a] and expand for X/ry) <1 


xX 2X 
+ *) + ro po)? = — +4, [5a] 
Td To 
3 Integrate Equation [3b] and expand for X/r, < 1 
X 
(14%) ade +H [5b] 
rb 


4 Assuming r, and a, to be constants, Equation [5b] 
may be integrated from 0 to 7 and substituted into [5a] 


2 
To To To 
B Cc 


Using the fact that X = 0 = X at r = 0, Equation [6] may 
be integrated to obtain 


Expanding this solution into a series but collecting only the 
second- and third-degree terms in +r (terms beyond 7° are 
affected by the assumptions made earlier), we have 

H 
2 r,? 3 8] 


Now, from Equations [4b and 5b] we have 


dg, H rar+H 
dr rp? dr re(1 + 2X/rs) 


By integrating the difference between Equations [9] 
52 
rd — Xar aD) 


Utilizing Equation [7] 
H 


— 
2 a 


C  cosV/Br 1 a, siny/ Br 
ve) 
Expanding this solution and collecting second- and _ third- 
degree terms 


— do) = 


[12] 


2 
3 


Thus, Equations [7 and 11] in trigonometric form and 
Equations [8 and 12] in series form present the radial and 
circumferential deviations from the ballistic reference orbit. 

The position of the vehicle has now been obtained; also 
useful are two additional parameters: The vehicle velocity V 
and the inclination angle y. See Fig. 2. 

Velocity may be obtained from an energy balance. The 
total energy at time 7 equals the energy of the ballistic ref- 
erence orbit plus the subsequent work done by the radial and 
circumferential components of thrust as follows 


y? 1 1 


Again neglecting higher degree terms, this equation becomes 


2 r 2 To i 


The parameter y may be determined from a consideration of 
angular momentum. The angular momentum at time r+ 
equals the angular momentum of the ballistic reference orbit 
plus the angular momentum produced by the circumferential! 
component of thrust as follows 


rV cosy = H+ a. rdr’ [15] 


Substituting from Equation [8] and neglecting higher degree 
terms we have 


y = cos"! (H + aersr)/rV [16] 


Summarizing, the parametersr, ¢, V and y have been deter- 
mined as functions of r. The parametric assumptions made 
in the derivation were that a, and a, = constant < 1, that 
X/r, <1, and that r, = constant. 
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Fig. 3 Approximate and exact characteristics of departure tra- 
jectories from Earth’s orbit for a, = 3 anda, = 0 


Comparisons of Approximate and Exact Trajectory 
Characteristics 


Analytic solutions to the general equations of motion (Eq. 
[3]) cannot be determined. However, Copeland (3, 4), has 
determined analytic solutions for the case of radial thrust 
only. His analytic solutions fora, = } (this value is approxi- 
mately 0.8 X 10-4 Earth surface g) will now be utilized for 
comparison with the foregoing approximate solutions. 
| Exact and approximate solutions of r as a function of ¢ are 
presented in Fig. 3. The dashed line approximations which 
range from 0 < t < 150 days are solutions of Equations [7 and 
8], where the reference ballistic orbit is taken as Earth’s orbit 
about the sun (rp = 1 forallt). It may beseen that att ~ 79 
days the differences between the approximations of Equations 
[7 and 8] and the exact solution are Ar ~ 0.002 and Ar = 
0.014, respectively. 

Another ballistic reference orbit, specified by the require- 
ment that its point of osculation coincide with the exact orbit 
atr = 1.10 (t ~ 79 days), was determined from the equations 
of (1) and is also presented in Fig. 3. Approximations from 
Equations [7 and 8] are presented in the range 79 <t <200 
days and were found to be more accurate than those corre- 
sponding to the previous range. For example, at t ~ 158 

days (time from osculation plus 79 days) the differences be- 
tween the approximations of Equations [7 and 8] and the 
exact solutions are only Ar = 0.001 and Ar = 0.005, 
respectively. 

Figs. 4 and 5 present an important application of the fore- 
going approximations. The fact that new reference orbits 
may be periodically determined from the approximate equa- 
tions permits the calculation of trajectories which closely 
agree with the exact solutions for long periods of time. This 
procedure was used to calculate a trajectory for the distance 
to Mars. 


‘In Equation [7] a constant value of r, was assumed. The 
value chosen was approximately an average value for the time 
interval of interest. 
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Experience gained from the calculations for Fig. 3 dictated 
the use of Equation [7] because of its initial accuracy for the 
departure phase of the trajectory and the use of Equation [8] 
because of its extreme simplicity for the remaining phases. 
New reference orbits were calculated at t ~ 79 days andt ~ 
153 days. The results, r and ¢, are presented as a function of 
tin Figs. 4 and 5, respectively. 

The deviations at the distance of Mars (r ~ 1.52) per pa- 
rameter interval for the time interval 0 S ¢ S 253 days were 


A, Ad 
152-100 0.8 per cent 1735-0 0.9 per cent 
From another viewpoint, the deviation in time between the 
approximate and exact solutions at r = 1.52 is approximately 
2.3 days out of 255.3 days. Although not presented, excel- 
lent comparisons were also obtained for V and y. 

These results are quite suitable for feasibility studies or 
preliminary engineering estimates. Although these results 
were obtained for the case of radial thrust only, it is felt that 
similar results could be obtained for the general case of con- 
stant thrust at an arbitrary angle 6. Unfortunately, no 
accurate numerical solutions of these latter cases were avail- 
able for comparative purposes. 

The equations of this paper can also be applied to other 


types of trajectories. For instance, trajectories utilizing vari- 
able thrust programs may be approximated by a sequence of 
steps of constant thrust at a given angle. Another type of 
trajectory for which the simple approximate methods are 
ideally suited is a ballistic trajectory which utilizes application 
of low thrust for corrective maneuvers. 

The foregoing results are not limited to heliocentric applica- 
tions. In fact, planetocentric maneuver or escape trajectories 
may be obtained from the previous calculations of heliocentric 
trajectories by merely changing the gravitational constants. 
For example, if we wish to consider geocentric orbits, we need 
only let Ro equal 1 Earth-radius and go equal the gravitatioual 
acceleration of the Earth at Ro. 

The authors wish to express appreciation to Dr. A. Kartvcli, 
Vice President of Republic Aviation Corp., for permission to 
publish this analysis. 
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Normal Dispersion of a Re-Entry Body 
R. A. NORLING! 


Bendix Corp., Advanced Development Lab., South Bend, 
Ind. 


A solution is obtained for the normal dispersion or dis- 
placement of a re-entry body from its static path due toa 
small aerodynamic trim misalignment for the planar case. 
An expression is obtained which relates the normal dis- 
persion as a variable function of altitude. The solution 
has been further simplified for selected parameters of the 
ballistic factor 8 and re-entry path angle yz. 


ONSIDER a nonspinning re-entry body on the straight 
line trajectory shown in Fig. 1. Furthermore, assume 
that the body possesses an inherent trim misalignment ar. 
Since it has been shown in (1)? that the effect of oscillatory 
terms on the normal displacement is negligible provided that 
the re-entry body is relatively stable, the angle of attack for 
our purposes can be considered constant a = az, i.e., equal 
to the trim angle. 

As the body penetrates into the sensible atmosphere, the 
straight line nature of the trajectory will be altered by virtue 
of the perturbing lift forces resulting from the small trim 
angle ar. The base or static velocity can be represented by 
the expression derived in (2) 


V= Vge~* 
a = static trajectory parameter® = ps.1.g/2cB8 sin yz [2] 


(Vz = re-entry velocity) {1} 


p = atmospheric density ratio, p/ps. 1. = e~™ [3] 
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For the planar case, the transverse equation of motion for 
small trim angles becomes 


—mVY = (1/2)pps. 1.V*CLgarA (4] 


In this equation, y is considered as the deviation from a 
straight line trajectory, and along with its derivatives, to be 
of perturbation order. Equation [4] can be rewritten in the 
form 


-VY = as pV? [5] 


Measuring time from the initiation of re-entry, the velocity 
normal to the trajectory and in the plane of the resulting lift 
force becomes 


i, ) 28. Jo? 


dp dpdy dp V si 
— =—— = ——Vsin 
dt dy dt dy id 


which when combined with Equation [3] give 
dt = dp/cpV sin vz [6] 


Performing the change of variable, the normal velocity be- 


comes 
Cla 9Ps.- L- 
d 
@) 2cB sin yz J 0 
assuming that jp > 0 as y > high. Upon integration the 
preceding equation yields 


= (KVe/a)(1 — e~%) [7a] 


where 


Cla Ps. L- 
K = 7 
2cB sin ve [7] 
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Similarly, the normal displacement becomes 


which upon integration and substitution of limits yields 


ap d(ap) _ 
Kab) = | (nap In [8] 


The first term which appears in the brackets may be re- 
written in the following manner 


0 


0 ap ap Fata ap 
+ © daz 
ap 


where Ex(ap) is the exponential integral of argument af. 
Sub-tituting Equation [9] into [8] and recognizing that the 
Eul:r constant may be written as (3) 


= lim + In = —0.5772 


z—0 


We obtain the result 


K 
&(ap) = ———— [Ei(ap) — In ap — ¥] 
ac sin Yr 


Substituting the appropriate values for K and a (Eqs. [2 and 
7b]), we obtain the expression 


t(ap) = [Ei(ap) — In ap — y] 


Cp c sin YE 


Similarly, we may rewrite the expression for the normal 
velocity (Eq. [7b]) as 


E(ap) = (1 ) (11) 
D 


Equation [10] expresses the normal deviation from a straight 
line trajectory for a small trim condition in terms of tabulated 
functions. The equation is subject to the restrictions im- 
posed by the assumption that re-entry occurs along a straight 
line, i.e., gravity is neglected (2). 
Simplified Approximation 
Equation [10] can be rewritten in terms of an infinite series 
by noting that (4) 
Ei(ap) = 0.5772 + In (ap) + ap + 
(ap)? (ap)* (ap)* 
2!2 3!3 4!4 


then 
Ei(ap) — In (op) — 0.5772 = ap | 1 + 4 ] 
2!2 3/3 


For ap < 0.4, the error is 11 per cent or less if all higher degree 
terms are neglected. Equation [10] then becomes 


Cie ar 
(<=) (. sin —) [12] 
assuming that 
= 0.0027 slug/ft? 
g = 32.2 fps? 
c = 1/23,500 ft- 


we obtain a = 1020/8 sin y. Therefore, for the approxima- 
tion [12] to be valid to within 11 per cent, the following in- 
equality must be satisfied 


B sin Yzx/p > 2540 [13] 
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Fig. 1 Re-entry geometry 
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Fig. 2 is a graph of the inequality [13] indicating the region 
where the approximate expression [12] is valid. 
In terms of the trajectory parameters, Equation [12] be- 
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An approximate expression can be similarly obtained for 
the velocity expression [11] by expansion of the exponential 
term 


b= (G2) arve[1 - (1-9 -+...)] 


Letting a = 1020/6 sin yz, the expression for the normal 
velocity becomes 


Cia arVe \_ 

= 1020 (=) (42 
Equations [14 and 16] represent simplified expressions for 
the normal displacement and velocity histories in the 
plane, which include the resultant lift and drag forces. 
These forces may not be coincident with the plane of the 
trajectory, and, therefore, the normal dispersion described 
here will, in general, manifest itself as a combined range and 
azimuth error to the overall missile trajectory. The dis- 
persion described here may be greatly attenuated by spinning 
the missile about its longitudinal axis, thereby averaging 
out the lift perturbation about that axis. For most bodies, 
the constant nature of the aerodynamic quantities C,, and 
Cp are maintained for continuum flow for Mach numbers 
in excess of 4. The altitudes at which a certain velocity or 
Mach number occurs remain a function of 8 by virtue of 
Equation [1]. 


Nomenclature 


static trajectory parameter, sin 
characteristic area, ft? 

atmospheric density term, 1/23,500, ft=! 
drag coefficient, assumed constant 

lift curve slope (1/radian), assumed constant 
acceleration of gravity 

defined in Equation [7a] 

mass of re-entry body, slugs 

time, sec 

initial re-entry velocity 

velocity, fps, Vze~*” 

weight of re-entry body, lb 

altitude, ft 

trim angle, radians 

ballistic factor, W/CpA, assumed constant 
re-entry path disturbance, Euler constant 
initial re-entry path angle 

displacement normal to the static re-entry path 
atmospheric density ratio, p/ps.x. 
atmospheric density, ps.z.e~°”. [See (5)] 
ps.t. = sea level density, 0.0027 slug/ft® 
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Structural Effects and Particle Content 
of Interplanetary Space 


L. REIFFEL! 


Armour Research Foundation, Chicago, III. 


Calculations on structural damage by solar plasmas are 
presented which show that solar helium is the dominant 
influence in sputtering, whereas solar hydrogen dominates 
radiation damage effects. Newly available laboratory data 
are used to calculate erosion rates in space. 


NOWLEDGE of the particle content of interplanetary 

space is important to our understanding of a wide 
variety of astrophysical phenomena. Information on densi- 
ties, velocity distributions and chemical composition of the 
interplanetary gas also assumes significance in evaluating the 
effect of the environment of space on manmade structures, 
especially those involving thin sections or surface films. We 
have treated some of these effects previously (1).2. Estimates 
on most important parameters of the interplanetary medium 
have been largely inferential, as discussed in (1). Generally,* 
protons and electrons are considered to be present in equal 
number densities ranging from 10? to 10% per cm’ under quiet 
solar conditions to perhaps 10° per cm# in intense streams 


Received Feb. 18, 1960. 
of Physics Research, Physics Division. Member 
? Numbers in parentheses indicate References at end of paper. 
3 We take this opportunity to note some pertinent references 
‘3 through 6), accumulated since the submission of (1) to ARS 
OURNAL. 
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issuing from an active sun. Corresponding velocities, sug- 
gested by time correlations between solar and terrestrial 
events as well as other data, appear to lie between 500 and 
1500 km per sec or more. The magnitude of the velocity 
under quiet conditions is unfortunately particularly uncertain. 
Note that a proton velocity of 1000 km per sec corresponds 
to an energy of 5 kev. 

Direct measurements on the density of particles in inter- 
planetary space have apparently not yet been carried out, 
with the exception of one set of observations made by the 
USSR. The Russian results indicated densities fluctuating 
between 10? and 10° particles per cm’. It is interesting to 
note that these fluctuations did not show any obvious cor- 
relation with interplanetary magnetic field fluctuations on 
which observations were also made. Attempts to measure 
the velocities of the particles in the interplanetary gas were 
not successful.‘ 

In our earlier work, both radiation damage and sputtering 
effects on thin films or surface coatings for average conditions 
and encounters with intense streams issuing from an active 
sun were considered. In treating sputtering effects, we 
originally employed an average sputtering yield 6 (8 = num- 
ber of atoms per incident ion) of unity, while explicitly noting 
the uncertainty of the value and the simplicity of scaling to 
other values of yield as more complete data became available. 
Recently, Grgnlund and Moore (7) have made an extensive 
study of the sputtering of silver by hydrogen and helium using 


‘ Information on the USSR measurements of particle densi- 
ties, magnetic field fluctuations and attempts at velocity me:s- 
urement was provided for the author by L. I. Sedov and V. I. 
Krasovskii during conversations in Chicago, Nov. 20, 1959. Un- 
fortunately, information on Soviet estimates of the average 
interplanetary gas density was not obtained. 
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magnetically analyzed ion beams. Their results, for normal 
incidence, indicate substantially lower yields than have been 
obtained heretofore by gas discharge methods, such as em- 
ployed by Wehner (8). The discrepancy apparently arises 
from the presence of heavy ions in the gas discharge which 
are avoided by use of analyzed beams. The gas discharge 
technique can also be modified to eliminate the effects of 
undesired ions and may be combined with spectroscopic 
detection methods to provide a very sensitive means of 
measuring yields near threshold (9). Wehner’s measure- 
ments indicate that threshold values of about 50 ev are ap- 
propriate for most materials; yields in the neighborhood of 
threshold are 8 = 10~4 to 10-5, but rise very rapidly to values 
of the order of 8 = 10~! for energies > 150 ev. 


Grénlund and Moore’s detailed study of Ag covered the 
rang: 2 to 12 kev, in which range 6 exhibits a broad and flat 
maximum for light bombarding ions. At 5 kev, 6 for H* is 
0.035, whereas for Het, 6 is 0.48. Thus normal incidence 
sputtering by helium is 14 times that due to hydrogen, a ratio 
which is roughly maintained over the entire energy range 
studied. It is important to note that with the high 6 for 
He~ as compared with H*, the sputtering of objects in space 
will actually be dominated by solar helium (10). This 
rather surprising result occurs if helium travels along with 
hydrogen in the outward streaming solar gas and if relative 
abundances characteristic of the solar atmosphere are re- 
tained in the streams. Aller (11) has recently suggested a 
relative atomic abundance for H/He of 7/1 as the best esti- 
mate available. 


The overall effect of using the new absolute values for 6 
including the effects of solar helium is to give an effective B 
to be used with proton density figures of 0.1 at normal inci- 
dence. To account for increased yields at other than normal 
incidence [for some materials B at 30 deg is 10 times 6 at 
90 deg (9) ], a good average value for the effective 8 of a 1000 
km per sec solar stream is probably about 0.2 for materials 
similar to silver. This is one fifth the yield assumed pre- 
viously and therefore reduces the sputtering rates given in (1) 
for any particular proton density by a factor of 5. 


Were we to include effects of the heavy element admixture 
in the extended solar corona, using the same assumptions 
as we applied to He* and H* computations, the sputtering 
rates would be increased somewhat. Unfortunately, the gaps 
in our knowledge of the details of heavy element concentra- 
tion, acceleration and sputtering efficiency are all very great, 
so that even rough estimates of the erosion are difficult to 
make. However, besides their work with H* and Het, 
Grgnlund and Moore have measured yields for a variety of 
other ions bombarding Ag surfaces at 5 kev and find for N+ 8 
is 4.0, for O+ B is 4.4 and for Net 8 is 5.5. Thus, the effi- 
ciency of sputtering for ions of mass 20 is roughly 160 times 
greater than for protons. The average atomic weight for the 
heavy elements in the solar atmosphere may be taken as 
roughly 32, and the relative abundance, by weight, for these 
elements appears to be about 2 per cent (1) so their numerical 
abundance is of order 10~%. If we simply assume that the 
heavy components behave as 5-kev neon ions, then it follows 
that the sputtering they cause, when compared with that pro- 
duced by solar protons and helium, will result in an increase in 
erosion rates not exceeding 10 per cent. Ignoring this con- 
tribution tends to compensate for somewhat reduced proton 
densities which would be inferred from those observations 
which depend upon electron densities and charge neutrality 
arguments. 


Using 8 = 0.2, the calculated time for complete removal 
of an optically opaque, 300-A thick, silver film under bom- 
bardment by a 1000 km per sec solar stream containing 600 
protons per cm’, with accompanying heavier ions, is readily 
shown to be 1.5 X 10’ sec or 6 months. It is of interest to 
compare this result with etching rates derived from meteorite 
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cosmic ray age information, since in such a calculation one 
can avoid assumptions about absolute sputtering yields and 
particle densities. Whipple and Fireman (12) have re- 
cently suggested an upper limit etching rate of 30 A per year 
for an iron meteorite in an orbit near that of Earth. The 
calculation was based on a cosmic ray age, for the particular 
meteorite fall (Sikhote-Alin) involved, of 5 X 108 years to- 
gether with a cosmic ray attenuation length of 72cm. Work 
by Fisher and Schaeffer (13) using argon, helium and neon 
age data®> has meanwhile shown that the correct age for the 
meteorite is nearer 1.7 X 108 years. This has the effect of 
raising the etching rate limit to about 100 A per year for an 
evenly exposed spherical body. For an oriented iron surface, 
the erosion rate limit would then be 400 A per year. The 
relative sputtering rate for Ag is two to three times higher 
than that for Fe if results for bombardments with 400-ev 
argon ions are typical (9). Thus, the erosion limit, from 
meteoritic data, for an oriented Ag surface would be 800 to 
1200 A per year as compared with the result calculated here 
of 600 A per year. Uncertainties in meteoritic etching rate 
information are unfortunately considerable for reasons dis- 
cussed elsewhere (10) and because of questions concerning 
cosmic ray constancy in time and attenuation length in iron 
(14). 

The radiation damage calculations given in (1) remain es- 
sentially unchanged by the present considerations, since 
the particle densities are of the same order as used previously. 
Thus charring of plastics or other organic materials to depths 
of a few thousand angstroms in short times is to be expected. 
For thin films, such effects could be quite significant. Changes 
in inorganic surface coatings, as a result of color-center and 
defect production and similar radiation induced phenomena, 
are also to be anticipated. In addition to being a function 
of particle density, the importance of any of these possi- 
bilities is directly dependent on the velocity and hence the 
range of the incident ions. If actual velocities are much 
lower than those used in the computations, the smaller depths 
in which radiation induced changes would then be produced 
might make them tolerable. However, calculations using 
only a considerably reduced particle density figure (e.g., 
10 to 100 particles per em’) would still suggest appreciable 
radiation effects with many materials. Final conclusions on 
both sputtering and radiation damage effects must await 
direct measurements on the density, velocity and chemical 
composition of the interplanetary gas as a function of solar 
activity. 
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Magnetic Effect in a T-Tube 


BOYD W. HARNED! 


General Electric Co., Evendale, Ohio 


Data observed with five return lead variations in a T-tube 
electromagnetic plasma accelerator have been analyzed 
for a possible correlation between energy conversion ef- 
ficiency and geometrical configuration. An attempt has 
been made to separate out the magnetic portion of the 
total accelerative process. 


FURTHER analysis of data (1)? on the conversion of 
electrical energy to kinetic energy of air contained in a 
T-tube is given. 

Brass electrodes separated by a 30-mm gap in a 39-mm ID 
Pyrex T-tube were connected to a 10-uf capacitor source 
through a triggered spark gap switch. In order to determine 
the nature of the effect of the magnetic field at the gaseous 
discharge of the ordinary backstrap configuration (2), several 
other return lead connections were made to provide a varia- 
tion of magnetic field. Air pressure before discharge was held 
at 250 uw of mercury, and capacitor energy in all cases was 500 
joules. 

If the T-tube is visualized in Cartesian coordinates with the 
electrodes along the x axis and the desired magnetic field 
which repels the discharge along the y axis, then the exhaust 
arm of the T-tube will extend along the z axis. Ordinarily, 
the backstrap is parallel to the electrode alignment in the 
zx direction, and the y component of the field at the arc is pro- 
duced because of the circular pattern of the magnetic field due 
to the discharge current passing through the strap. Five 
cases will be reported here: 

1 The return lead located so as to minimize the magnetic 
field at the discharge—‘“no field’’ case. 

2 A length of }-in. copper tubing used as a backstrap in 
the zx direction. 

3 A two-turn coil wound of copper tubing, one on either 
side of the discharge section of the tube, so that the magnetic 
field is in the y direction. 

4 A four-turn coil similar to that of case 3, with two turns 
on either side of the tube and the magnetic field in the y direc- 
tion. 

5 Atwo-turn coil similar to that of case 3 but oriented with 
one turn in front and one behind the discharge in order to 
provide a field in the z direction and thus accelerate the 
gaseous conductor in the y direction—“wrong field” case. 

The data reported were: Frequency of the discharge ringing 
circuit, inductance and resistance calculated from the damped 
wave form observed with an oscilloscope, velocity of the ac- 
celerated gases as observed by a photomultiplier through an 
aperture system and a “Jight pipe,” and momenta of the gases 
as observed by means of a lightweight plastic disk pendulum. 
From these data were calculated an energy conversion effi- 
ciency value. The T-tube was shorted out completely so that 
the inductance and resistance of the external circuit could be 
observed and subtracted from the values in the various return 
lead cases, thus providing an estimate of how the added circuit 
inductance affected the action. Reproducible velocity data 
were obtained which showed a variation with distance along 
the exhaust tube, as reported by Kolb (2) and Kash et al. (3). 
All measurements referred to here apply to a distance of 18 
em from the discharge, accounting in part for the low ef- 
ficiencies. 

Shown in Table 1 are data for the five cases. The resistance 
values are exclusive of external circuit; the velocity values 
are those of the first shock only, although subsequent shocks 
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Table 1 Operating characteristics 


Configurations 
1 2 3 4 5 


Frequency, kc/sec 52.6 64.5 51.2 37.0 51.2 

T-tube discharge 0.024 0.032 0.064 0.083 0.042 
resistance, ohms 

Efficiency, percent 3.7 7.4 5.8 4.4 1.0 

Velocity, cm per 1.51 1.80 1.65 1.26 0.98 
microsec 

Energy dissipated 56 85 131 125 90 
in T-tube dur- 
ing first half cy- 


cle, joules 
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Fig. 1 Energy per unit resistance vs. efficiency for various 
T-tube configurations 


due to successive current maxima in the ringing contributed 
to the observed momenta. Dissipated energy values are cal- 
culated from observed circuit parameters. Resistance and 
inductance values of the gaseous discharge must be considered 
as average values throughout the energy dissipation process. 
T-tube resistance includes that of the electrodes, excluding the 
return lead. 

The equation of motion for a partially or fully ionized gas 
(4) subjected to combined electric and magnetic fields con- 
tains two terms of interest here, namely, a pressure gradient 
term and a j x B term. It seemed possible that the data 
would provide a means of separating these two terms in crude 
fashion, if a few simple assumptions might be imposed. Since 
ringing frequencies differed in all these cases, and hence cur- 
rent values and magnetic fields varied with time, an attempt 
was made to correlate the data in terms of energy dissipated 
in the T-tube during the first half cycle of the discharge. 
This seems reasonable, since it is assumed that the first shock, 
whose velocity is reported, is caused by the first half cycle of 
the current. Case 1, the “no field” case, served as a criterion 
for the pressure gradient term, by which the additional T-tube 
circuit resistances and dissipated energies observed for the 
other cases with magnetic fields were obtained by subtraction. 
A ratio was made of these energy differences divided by the 
resistance differences and was plotted against the efficiencies 
observed. Since the “wrong field” case 5 deflected the arc per- 
pendicularly to the exhaust direction and resulted in a reduced 
efficiency, this case was plotted as a negative effect. The 
resulting curve is shown in Fig. 1. 


(Continued on page 7085) 
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